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DEFINITION OF SYMBOLS AND UNITS 


Symbol 

Definition 

Units 

9 

attitude angle 

rad 

■8 

angle between reference and inertial 
velocity vector 

rad 

a 

angle of attack 

rad 

P 

control deflection angle 

rad 

V R 

velocity of vehicle relative to air 

m/ s 

V 

velocity of vehicle 

m/ sec 

Vw 

wind velocity 

m/ sec 

z 

direction normal to reference 

m 

z 

direction normal to vehicle centerline 

m 

X 

vehicle longitudinal axis 

m 

F 

total thrust of the vehicle booster 


m 

total mass of the vehicle 

kg-sec 2 /m 

I 

vehicle moment of inertia about the CG 

kg-m-sec 2 

X E 

engine moment of inertia about hinge 
point 

kg-m-sec 2 

D 

drag force 

kg 

S E 

engine first moment of inertia about 
hinge point 

kg-sec 2 

q 

dynamic pressure 

kg/m 2 

R* 

thrust normal slope of control engines 

kg/rad 

n ! 

aerodynamic normal force slope 

kg/rad 

C 1 

aerodynamic moment coefficient 

1/ sec 2 

n ' a 

aerodynamic bending moment coefficient 

kp-m/rad 

% 

control force coefficient 

V 

kp-m/ rad 



DEFINITION OF SYMBOLS AND UNITS (Cont'd) 


Symbol 

Definition 

Units 

c 2 

control moment coefficient 

1/sec 2 

a o 

displacement gain 

-- 

a l 

rate gain 

sec 

b o 

angle of attack gain 

-- 

§2 

accelerometer ga in 

sec 2 /m 

^a ~ x cg ” x a 

distance from vehicle CG to accel- 
erometer 

m 

£ cp = x cg ■ x cp 

distance from vehicle CG to the CP 

m 

£ cg = x cg - x h 

distance from engine gimbal to vehicle 

CG 

m 

h = x h - X E 

distance from engine gimbal to engine 
mass CG 

m 

m sk 

slosh mass 

kg-sec 2 /ml 

z sk 

slosh mass displacement, normal 
to reference 

m 

£sk 

slosh damping 

-- 

W sk 

slosh frequency 

rad/ sec 

x sk - *h 

distance from engine gimbal to 
slosh mass CG 

m 

-^sk _ x cg " x sk 

distance from vehicle CG to slosh mass 
CG 

m 

S 

cross sec tiona 1 reference area 

«,2 

m 

T ii 

generalized displacement of the i tb 
mode (usually denoted as "normal 
coordinates") 

m 

Mi 

generalized mass 

kg-sec 2 / m 
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DEFINITION OF SYMBOLS AND UNITS (cont'd) 


Symbol 

Definition 

Units 

Gj i 

bending frequency 

rad/ sec 

L 

bending mode damping 

-- 

Y i (x) 

normalized* displacement at Sta x 

-- 

Y^x) 

normalized* slope = [ Y - (x ) ] 

dx i- 

1 /m 

Y £ (x) 1 ( i 

displacement at Sta x due to i t ^ bending 
mode 

m 

Y ’i(x)\ 

angular displacement at Sta x due to i^ 1 
bending mode 

rad 

n< x >h 

angular rate at Sta x due to i^ 1 
bending mode 

rad/ sec 

Y i (x>*; t 

angular acceleration at Sta x due to 
ith bending mode 

rad/ sec 2 

Qi 

generalized force for the i^ mode 


W (x , t) 

force distribution over the length of 
the vehicle for all forces acting upon 
the vehicle 

kg/m 

m k 

mass 

kg - s e c / m 


moment of inertia 

kg-sec^ 

\ “ X h 

distance to the vehicle Sta x^ as 
measured from the vehicle gimbaL point 

m 

J 'k 

distance between two vehicle stations 

m 


damping 

-- 

L \ 

angular frequency 

rad/ sec 


* The bending displacement 
gimbal point. The slope 
men t . 


is normalized to "+1" at the vehicle 
is a function of the normalized displace- 
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TECHNICAL MEMORANDUM X-53036 


CONTROL THEORY HANDBOOK 
SUMMARY 


This report presents under one cover, the equations of motion and 
the basic control theory applicable to stability and response analyses 
for a flexible launch vehicle, using a unified coordinate system and 
notation. To provide some background and insight in the control of 
large flexible boosters moving through the earth’s atmosphere, five of 
the basic control problems are discussed. 

The flight system coordinates and notations are shown and the rigid 
body equations are derived for both the pitch and yaw planes. A con- 
ventional control system is introduced which containsa position gyro, 
rate gyro, accelerometer and angle of attack meter. The gains of the 
control mechanism and the vehicle parameters are related to the frequency 
and damping of the rigid body. Both the "Drift Minimum" and "Load Mini- 
mum" control principles are developed. 

The bending and slosh equations are derived by writing the energy 
expressions and then applying Lagrange’s equation. The method of com- 
puting bending modes and frequencies for a flexible body is shown for 
both a simplified continuous mass model and a lumped mass model. 

The construction of a synthetic wind profile for control system 
studies using the 95 or 99 percent probability of occurrence wind speed 
profile and the 99 percent probability of occurrence wind shear envelope 
is discussed and illustrated, and the method for superimposing a gust on 
the synthetic wind profile is shown. 

Block diagrams and the Laplace transform are introduced to relate 
the system equations in a form which can be studied in terms of general 
feedback theory. 

Routh’s stability criterion, Hurwitz’s stability criterion, root 
locus, frequency response methods, and Nyquist criterion are discussed 
and are applied to a vehicle containing one bending mode, a control 
filter and an actuator. The corresponding root locus plot, Bode plot, 
Nyquist plot and Nichols plot are drawn. 

The basic elements of an example adaptive control system are dis- 
cussed and its corresponding block diagram shown. 



The appendices contain the block diagram and transfer functions for 
several sensors and engine actuator. A summary of the flexible body 
equations, including the effects of engine inertia, bending motion and 
slosh motion, and a derivation of the bending moment at any station 
along the vehicle longitudinal axis are also given in the appendix. 


I. INTRODUCTION 

During initial design phases of a vehicle, numerous trajectories 
are computed using various degrees of sophistication. Almost invariably 
it is assumed that the vehicle is controlled perfectly and does not 
deviate from the desired trajectory. Once the trajectory has been 
defined, it is then the job of the guidance and control groups to 
specify a means of controlling the vehicle response in order to mini- 
mize errors in position and velocity. 

The overall guidance function is to determine the position and 
velocity of the vehicle and reduce these variables to a pitch and yaw 
command. These pitch and yaw- commands are then the reference inputs to 
the control system. 

The basic control function is to produce the desired output of a 
variable based on the reference input to the control system. This can 
best be illustrated by the functional block diagram shown below. 



FIGURE 1. FUNCTIONAL BLOCK DIAGRAM OF CONTROL SYSTEM 
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Here the reference input to the control system is R, and the 
desired output is C. The system contains the necessary sensors to 
detect the desired output and may be in the form of attitude, attitude 
rates, etc. The comparison between R and the feedback signal B results 
in an actuating signal E that is the difference between these two 
quantities. The actuating signal produces an engine deflection 0, 
which correspondingly produces the desired output. This system is 
called a closed-loop control system since it compares the output and 
input quantities to maintain the output at a desired value. 

The control system must contain the necessary sensors to detect 
the translational, rotational and vibrational motion of the vehicle. 

The commonly employed sensors are the rate gyro, position gyro, acceler- 
ometer and angle-of-at tack meter. 

The translational and rotational motion are conveniently controlled 
by deflecting the thrust vector or by use of movable surfaces on aero- 
dynamic fins. When aerodynamic control surfaces are used, jet vanes 
must be included in the rocket engine exhaust to produce the control 
force until the dynamic pressure builds up enough for the aerodynamic 
control surfaces to become effective. Thrust vector control is more 
desirable than fins due to a weight advantage. The vibrational motion 
consists primarily of fuel slosh and structural bending. This motion 
can be easily sensed, but is difficult to control. The vibrational 
nature of the vehicle is specified in the early design phases. The 
structural bending characteristics are determined by the structural 
stiffness and xsreight distribution and cannot be simply altered once the 
vehicle is manufactured. The fuel slosh characteristics are largely con- 
trolled by the tank geometry, internal tank baffles, and location of 
propellant tank relative to the center of gravity of the vehicle. 

In designing a control system for a flexible booster moving through 
the earth* s atmosphere, there are generally five major problem areas that 
must be considered. 


1 . 

Drift from reference 

trajectory. 

2. 

Aerodynamic loads. 


3. 

Structural bending feedback. 

4. 

Fuel slosh. 


5. 

Adequate response to 

d is turbances . 


The author acknowledges Mrs. Joyce Harmon for her contributions to 
the section on Methods of Stability Analysis and Mr. John Livingston for 
his contribution of the Section on Wind Representation. 
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II. DISCUSSION OF BASIC CONTROL PROBLEMS 


A. Drift From Reference Trajectory 

The drift of a vehicle from a desired trajectory can be detri- 
mental to the overall success of the mission. Most missions depend on 
the accuracy with which a vehicle can place the payload into a given 
volume of space with a predetermined velocity and direction. Since 
drift during the boost launch phase results in large position errors 
at later stages of flight, it must be held to an absolute minimum. 

The primary disturbance affecting drift is due to winds aloft. 
Other items which cause drift are thrust misalignments, gyro drift, 
center of gravity variations, etc., but these are of secondary impor- 
tance compared to w ind disturbances. 

Wind speeds can be of slowly varying nature, sudden gust or a 
combination of both. For a particular launch site, wind speed design 
profiles are determined based on numerous wind speed versus altitude 
measurements. Generally, 95 or 99 percent wind profiles are used in 
control system design and specify the maximum v/ind speed versus alti- 
tude for a 95 or 99 percent probability of occurrence. Also embedded 
gusts of varying wave lengths, peak-to-peak amplitudes and number of 
successive gusts can be superimposed on the wind speed profile to study 
the response of the control system. The wind speed profiles are con- 
structed by building up to the 95 percent or 99 percent maximum value 
using the 99 percent probability of occurrence wind shears as shown in 
Section VII. 


B. Structural Loads 


The aerodynamic forces imposed on a vehicle during its launch 
trajectory affect the dynamic response of the vehicle to control system 
commands. These forces may be significant from time of liftoff until 
the vehicle stages or emerges from the effective atmosphere. The magni- 
tude of these forces are dependent on the angle of attack, dynamic pres- 
sure, Mach number and the aerodynamic characteristics of the vehicle. 

Although the dynamic pressure (q) is not large during the 
initial phase of flight, the pitch program is usually initiated during 
this phase resulting in relatively large angles of attack (;), and, hence 
appreciable values aq. These loads caused by pitch over command can be 
critical when a step change in commanded attitude is programmed. The 
loads can be greatly reduced by changing the programmed command to a 
ramp function. The aerodynamic loading is typically most severe in the 
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high dynamic pressure region and occurs in the range of altitudes inhere 
winds and gusts are at maximum intensity. It is often necessary in the 
high dynamic pressure region to constrain the angle of attack to small 
values so that structural design limits are not exceeded. 

The magnitude of the aerodynamic force is obtained by integrat- 
ing the distributed loads due to pressure and viscosity along the length 
of the vehicle. The resultant aerodynamic force is then resolved into 
a normal force component perpendicular to the longitudinal axis and a 
drag component along this axis. The normal force component acts at the 
center of pressure and is one of the primary influences affecting the 
translational and rotational motion of the vehicle. 

The location of the center of pressure (c.p.) relative to the 
vehicle center of gravity (c.g.) is an important consideration as it 
places a requirement on the amount of control torque necessary to main- 
tain stability and control of the vehicle. The location of the c.p. is 
primarily a function of (1) the position of the cones, flares and fins 
which make up the vehicle configuration and (2) the Mach number. 

Vehicles having the c.p. forward of the c.g. are aerodynamica 1 ly 
unstable and tend to deviate from the direction of the relative velocity 
vector. In the case where limited control torque is available, it may 
be necessary to introduce fins at the base of the vehicle to shift the 
c.p. aft. The vehicle is aerodynamica lly stable when the c.p. is aft 
of the c.g. and tends to align its longitudinal axis along the relative 
velocity vector. Generally, fins are not required for this case. 

In the case where loads become a problem, it may be necessary 
to include a load relief feature in the control system. Some load 
relief systems are designed to operate only in the region where the 
lateral acceleration or the product of angle of attack and dynamic 
pressure, aq , exceeds some preset limit. The preset limit is dependent 
on the structural limitations of the vehicle. Once the preset limit 
is exceeded, gains are produced which are proportional to the lateral 
acceleration or angle of attack depending on the type of sensors utilized 
This gain then becomes a dominant part of the engine command signal. 

The effect of this added engine command signal is to turn the 
vehicle into the wind, thereby reducing the angle of attack and ulti- 
mately the engine deflection itself. The decrease of both of these 
variables reduces the bending moments which are directly proportional 
to the angle of attack and the engine deflection. 
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The most common load relief system considered in the past is 
designed to blend load minimum control gains in and out as a function of 
time. The time interval chosen generally corresponds to the maximum 
dynamic pressure region where aerodynamic loading is the most severe. 

Since load relief type control systems attempt to align the 
vehicle attitude along the relative velocity vector, there is no pre- 
dominant control of the vehicle attitude \ . Therefore, y can assume 
large values, and the vehicle will drift away from the desired trajectory. 
Application of load relief type control is acceptable for short durations 
provided a deviation from the desired trajectory is permissible. 

The bending loads imposed on the vehicle structure are depen- 
dent on the type of control law used since for each type of control law 
a different response in engine gimbal angle -3, and hence eg results. 
Appendix G shows that the bending moment for a rigid vehicle can be 
expressed as 


bending moment = (M^)g 4* (Mp)[3. 

The form of the ;3 response is dependent on the type of feed- 
back included in the control law. For an attitude control law, the 
engine response contains signals from both a position gyro (y) and rate 
gyro (g) and is of the form :3 = a Q y 4 a rg where a D is the attitude gain 
and a 1 is the rate gain. If angle of attack ( /) is included in the con- 
trol law, the p response is a linear combination of g, and (X is of 
the form 


P = a f 


+ a x cp + by 


where b Q is the gain on angle of attack. The method for determining the 
gains is discussed in Section IVD. 

Figure la [10] shows a comparison of the effect of three control 
laws on bending moments due to a triangular-shaped wind profile of unit 
amplitude. The maximum bending moment is plotted in each case versus 
the wind spike duration for a typically large vehicle at the time of 
maximum dynamic pressure. The wind spikes of short duration in this 
figure correspond to a high shear, high fundamental frequency wind. 

The longer duration wind spikes correspond to a low shear, low funda- 
mental frequency wind. 
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FIGURE la. EFFECT OF CONTROL LAW AND WIND SPIKE DURATION ON BENDING 
LOADS AT MAXIMUM DYNAMIC PRESSURE FOR A SATURN CLASS 
VEHICLE 


It is seen that attitude control produces lower bending 
moments for short duration wind spikes while both types of control 
employing a feedback are most effective in reducing Loads due to long 
gradual wind spikes. 


C. Structural Bending Feedback 

One of the problems associated with the control of a flexible 
missile is the control of the rigid body such that structural bending 
is not reinforced. This problem arises because of the fact that struc- 
tural bending modes may be excited by control actions required for 
maneuvering, or by aerodynamic loading due to sudden wind gusts. The 
elastic vibrations are sensed by the attitude, attitude rate, acceler- 
ometer and angle of attack sensors and are fed back into the control 
sys tem. 


Since these sensors are mounted on the vehicle elastic axis, 
they detect the resultant motion of both tlie rigid body and the elastic 
vibrations. The elastic vibrations comprise various frequency components 
and would ordinarily be included in signals to command engine deflection. 
Without some form of compensation of the elastic vibrations, there is 
the possibility that resonance would occur between the bending mode fre- 
quency components and the control frequency producing an undesirable 
amplification of the bending mode. 

Conventional methods of compensating for structural feedback 

are 


1. Placing an electrical filter network in the feedback 
loop which passes the portion of the sensor signal 

in the low band of control frequencies and suppresses 
the higher frequency signal which can be associated 
with the elastic vibrations. 

2. Location of feedback sensors to minimize excitation 
due to structural vibrations, e.g., mounting the 
rate gyro near an antinode of the first bending mode, 
or mounting the accelerometer at the first mode node. 

3. Placing notch filters in the feedback loop to suppress 
sensor signals in a selected small range of frequencies, 
the center of the range or "notch” being the frequency 
of the bending mode. 
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Another item of importance in considering the structural bend- 
ing of a vehicle is obtaining accurate mode shapes and slopes of the 
bending modes. The complex structure (multiple tanks and engine clusters) 
and large size of some boosters make it extremely difficult to obtain 
accurate theoretical and experimental data describing the vehicle’s 
elastic properties. The modes shapes are therefore inaccurately known 
and vary considerably during flight because of mass distribution changes 
from propellant depletion. 


D. Fuel Slosh 


Some of the basic parameters which influence the stability and 
control of a vehicle are (1) the magnitude of the slosh mass*, (2) the 
frequency of the slosh mode, and (3) the position of the slosh mass along 
the longitudinal axis. Since more than 90 percent of the total weight of 
the vehicle is liquid, special attention must be given to the influence 
of fuel motion on the stability of the vehicle. 

The magnitude and frequency of the slosh mass are largely a 
function of the slenderness ratio >i/d (long th- to-d iame ter ratio) of the 
tank. For long cylindrical tanks, i.e., large slenderness ratio, the 
ratio of the slosh mass to total mass of fluid in the container is small, 
while for short tanks a large percentage of the mass is sloshing. 

For cylindrical tanks with large slenderness ratio, the slosh 
mass for the first mode remains constant as the propellants are depleted 
from the tank until the fluid surface is about one diameter in height 
from the bottom of the tank. Thereafter, the slosh mass decreases 
rapidly. The magnitude of second mode slosh mass is less than 3 per- 
cent of the first mode slosh mass and the mass of successively higher 
slosh modes is even less. In most practical cases, therefore, all 
slosh modes higher than the first can be neglected. 

Tanks having large diameters, i.e,, those with small slender- 
ness ratio, have two undesirable characteristics: the slosh mass 

increases and the natural frequency decreases. The increased slosh mass 
magnifies instability depending on its location along the vehicle longi- 
tudinal axis. Considering only a rigid vehicle with ideal control and 
one propellant tank, a danger zone for instability occurs when the 
slosh mass is loacted between the c.g. of the vehicle and the instantaneous 


"A portion of the total liquid mass assumed to move as an equivalent 
spring-mass-damper system. The magnitude' of the slosh mass is dependent 
on tank geometry and fluid density. 
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center of rotation.* The instantaneous center is on the order of 10 
to 15 meters forward of the c.g. for a typically large vehicle. For 
increasing slosh mass in this region, more damping is required to 
guarantee stability and the maximum amount of damping is required when 
located at the midpoint between the instantaneous center and c.g. This 
is generally the region where maximum bending displacements occur also. 
Considering the elastic vehicle, the danger zone for instability 
increases both forward of the instantaneous center and aft of the c.g., 
and also more damping is required as compared to the rigid vehicle. 
Ideally, the stability can be improved by locating the slosh mass in 
the region aft of the c.g. if there is a large separation ratio between 
the slosh frequency and the control frequency. 

Tanks with large diameters exhibit low natural frequencies. 

This tends to restrict the choice of the control frequency since the 
slosh frequency should be as far above the control frequency as possible 
to prevent resonance. As would be expected, the forces and moments on 
the tank due to fluid motion are magnified considerably near the resonant 
frequency of the fluid. Subdividing a large tank into smaller containers 
may be necessary to increase the separation between the slosh frequency 
and control frequency. This also reduces the slosh mass, which is a 
more important effect. 

The addition of baffles in the tank is effective in increasing 
the damping of the fluid motion. This method is commonly used to improve 
stability even though there is an associated weight penalty. 

Control sensors can enhance the stability of the vehicle con- 
siderably. Proper selection of type, location, gain values and vibra- 
tional characteristics can be used to suppress the hazard due to fuel 
slosh. 


E . Adequate Response to Disturbances 

A primary requirement of the control system in addition to 
maintaining stability and control of the vehicle is to maintain a desired 
transient and steady state response for any given disturbance. The pri- 
mary disturbance is due to winds aloft and may vary randomly in speed, 
direction, and frequency. 


Instantaneous center of rotation - as used in the discussion above the 
instantaneous center of rotation is taken with respect to the engine 
gimbal point and is the point (x^ c = x C g + Z/ cp) on the longitudinal 
axis about which the vehicle has only an angular velocity. 
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The control system must be designed such that the transient 
response to a disturbance does not impose excessive loads on the vehicle 
as would be the case for a very fast responding system, and yet not so 
sluggish that large errors accumulate in the vehicle 1 s attitude, angle 
of attack, engine deflection, etc. Two parameters which are indicative 
of the transient response are the natural frequency, oj n , and damping 
ratio, The values of :,; n and £ for the rigid body transient generally 

are chosen, based on peak overshoot and settling time requirements. The 
time of peak ove rshoot occurs when the product u^t = it radians, where 
J d = ‘h is the damped frequency. The maximum peak overshoot 

for a second order system due to a step input is 


Maximum Peak Overshoot = lOOe 



percent 


and is only a function of £. 'In most control systems, a peak overshoot 
of 5 percent is desirable and corresponds to £ = 0.707. Accordingly, 
the damping ratio £ may be called a figure of merit, which identifies 
the maximum overshoot of the system variables due to wind inputs. 

The settling time for the transient response of the rigid body 
mode can be identified with its natural frequency. The system is con- 
sidered to have reached steady state when the response remains within a 
specified tolerance band. The settling time is defined as t s = k/to n , 
where k is obtained from the transient response curve as shown in Fig- 
ure 2. The value of k will vary depending on the value of damping 



FIGURE 2. IDENTIFICATION OF PEAK OVERSHOOT AND SETTLING TIME 
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ratio (£) or natural frequency (oj n ) chosen. For a given iu n , the settling 
time becomes longer as f is decreased due to the larger overshoots in 
the response. The effect of on settling time can best be seen by 
comparing two different systems having the same f. If the response 
shown in Figure 2, which is plotted in nondimens ional coordinates, 
applies to two different systems having the same then the system 
having the largest natural frequency will have the shortest settling 
time . 


The natural frequency of the rigid body mode should be kept 
sufficiently low so that there is a wide separation with frequencies of 
other degrees of freedom. Generally, the first slosh mode frequency 
establishes the upper limit on the control frequency of the rigid body 
mode since the first slosh mode contains the lowest frequency component 
as compared to the other degrees of freedom. The control frequency for 
the rigid body mode should be selected such that it is sufficiently 
below the first slosh mode frequency so that strong coupling does not 
exist and yet not so low that long settling times result. 

The transient and steady state response curves for a large 
vehicle are shown in Figures 3 through 6 for a step input wind distur- 
bance of 44 m/s. This corresponds to a wind angle ( > v of 5° as shown 
in the illustrations. These figures illustrate the effect of damping 
on the overshoot and settling time for this particular vehicle. The 
control frequency co n is held constant for each case of .30 cps and the 
damping ratio is varied from .30 to .80. 

The curves in Figures 3 and 4 indicate that the vehicle with the 
damping ratio ranging from .3 to .5 contains oscillations which persist 
for a relatively long period of time. For this low damping case, the 
overshoots are generally large and the settling time of the transients 
are around 6 seconds. Increasing the damping produces two desirable 
effects: it decreases the overshoot and reduces the settling time. 

The curves in Figures 5 and 6 show a more desirable response 
since they contain less overshoots and shorter settling times. The 
transients in these figures make one overshoot and then settle out to 
the final steady state values. The steady state values in each of the 
above four cases are seen to be independent of damping for this particular 
system with the exception of z, which continues to increase with damping. 

Increasing the C^/ C ratio (aerodynamic to control moment coef- 
ficient; see Section IV) has the effect of increasing the overshoot on 
the transient response. This is primarily due to an increase in the 

magnitude A for a sinusoidal transient of the form Ae J ^ sin (uj^t + ; ). 
Smaller vehicles generally have a large C ± / C ratio which decreases as 
the size of the vehicle increases. 
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FIGURE 3. TRANSIENT RESPONSE CURVES FOR A RIGID BODY 
DAMPING OF .3 AND FREQUENCY OF .3 CPS 



FIGURE 4. TRANSIENT RESPONSE CURVES FOR A RIGID BODY 
DAMPING OF .5 AND FREQUENCY OF .3 CPS 



FIGURE 5. TRANSIENT RESPONSE CURVES FOR A RIGID BODY 
DAMPING OF .7 AND FREQUENCY OF .3 CPS 
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FIGURE 6. TRANSIENT RESPONSE CURVES FOR A RIGID BODY 
DAMPING OF .8 AND FREQUENCY OF .3 CPS 



The natural frequency of the rigid body mode should be widely 
separated with frequencies of other degrees of freedom. The guidance 
mode establishes the lowest frequency existing in the system. The 
period of this mode is several orders of magnitude longer than the con- 
trol and elastic modes. The control mode in this case can be considered 
as uncoupled from the guidance mode. The upper bound on control fre- 
quency is established by the first slosh mode frequency or the first 
bending mode frequency whichever is less. The control frequency for 
the rigid body should be chosen so that it is sufficiently below the 
upper bound frequency to reduce mode coupling effects, and yet not so 
low that long settling times result. 


III. FLIGHT SYSTEM COORDINATES 


The orientation of the missile axes in inertial space is shown in 
Figure 7. 


Y 



The XYZ- coordina te system is defined relative to the reference tra- 
jectory as follows: The X axis is directed along the desired velocity 

vector V, the Z axis is normal to the plane of the desired trajectory, 
and the Y axis is in right-hand relation to X and Z. The vehicle orien- 
tation is defined by the xyz-coordina tes where the x axis is directed 
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along the vehicle longitudinal axis, the y axis is directed along the 
vehicle yaw axis, and the z axis is directed along the vehicle pitch 
axis. The origin of the xyz-coordinate system is located at the vehicle 

c.g. 

The pitch attitude Op and yaw attitude q,y define the direction of 
the missile longitudinal axis in the XYZ-coordinates . These angles are 
assumed small so that the equations of motion may be linearized. 

The degree of freedom along the X axis is eliminated by allowing 
the coordinate system to accelerate with that of the vehicle c.g. in 
the X direction. The equations of motion, therefore, will allow only 
accelerations relative to the XYZ-coordinate system in the Y and Z direc- 
tions. The orientation of the coordinate system is chosen in this man- 
ner since we are only interested in the attitude deviations and position 
and velocity deviations from the reference trajectory. 

Figure 8 below shows the orientation of the vehicle relative to 
the reference trajectory and the inertial coordinates XY. This is the 
geometry used in making a fixed time point study for a control system 
in the pitch plane. In this case, it is permissible to reduce the 
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problem to two dimensions since the cross coupling between the yaw plane 
and pitch plane is assumed small. Also shown in this figure are the 
angles X and X c , where X c is the attitude commanded by the pitch program 
and X is the actual attitude of the vehicle. The angle X c is obtained 
from the reference trajectory and is input to the vehicle autopilot. As 
can be seen from Figure 8, the pitch attitude * = X - X. The variable 

Op used in the control system equations rather than X c and X since it 
represents a small angular deviation from the desired directions X. 

Figure 9 below shows the orientation of the vehicle relative to 
the yaw plane of the reference trajectory. 



i 

Plane of 
Reference 
Trajectory 

4 


►- 

f 



m 

K 


t 


FIGURE 9. YAW PLANE GEOMETRY 


Again the problem is reduced to two- dimens ional motion in the yaw 
plane by assuming the cross coupling between the pitch plane and yaw 
plane is small. 
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IV. RIGID BODY EQUATIONS 


A. Derivation of Equations 


The free body diagram for a rigid vehicle moving through the 
earth's atmosphere is shown in both the pitch and yaw planes of flight 
in Figure 10 and Figure 11, respectively. The equations of motion are 
similar in both the pitch and yaw planes of flight; the discussion which 
follows will be for the pitch plane of flight. Subsequently, it will be 
shown how the pitch plane equations can be reduced to represent the yaw 
plane equations. 

The equations of motion are written for a nonrotating earth. 

In this case V represents the instantaneous velocity of the vehicle rela- 
tive to earth. The vehicle attitude is commanded by X c to align along 
the X direction to maintain the desired velocity . Due to the wind 
disturbance V w , the attitude error is cp = X c - X and the resulting 
velocity error is given by V G . 

The angle 7 denotes the bias of the reference direction X 
(direction of the nominal acceleration) from the desired velocity . 

For a gravity tilt trajectory the commanded rate X c is adjusted to cancel 
the component of gravita tiona 1 acceleration normal to V. In this case 
X lies along and 7 is zero. 

The coordinate system (X, Y) is oriented relative to the local 
vertical at launch by the commanded attitude X c . The vehicle x axis is 
oriented relative to the X coordinate by the attitude error cp. All 
forces acting on the vehicle are positive in the x and y directions. 
Moments about the vehicle center of gravity are positive when they pro- 
duce angular accelerations in the positive direction' of cp. The angles 
cp, a, a v , 'd and 7 are assumed small and the angle ■}( is negligible during 
the entire boost phase of flight. All angles are positive in the pitch 
plane for a counterclockwise rotation. 


The angle X c is used with the wind velocity V w and flight 
velocity V to determine the wind angle of attack o^. When computing the 
wind angle, the assumption is made that 7 = \|r = j = 0 and V w « V. The 
geometry for the condition is shown in Figure 12. The wind angle is 
defined as follows: 


a 


w 


tan - * 1 



cos X 


V w sin X c 


( 1 ) 


In the yaw plane the angle X c is zero and the wind angle for 
Vw « V becomes 


a 


w 


tan ” 1 


V 

w 

V 


V 

w 

V * 


1 n 
I o 


( 2 ) 




FIGURE 10. RIGID BODY COORDINATES, PITCH PLANE 





FIGURE 11. RIGID BODY COORDINATES, YAW PLANE 
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Local Vertical 



FIGURE 12. WIND ANGLE OF ATTACK 0^ 


The acceleration of the vehicle is obtained by finding the 
rate of change of velocity in the (X, Y) coordinate system. The 
expression for the velocity of the vehicle is 


V = V cos -Qi + V sin 


(3) 


The acceleration then becomes 

A = — = V cos di - V sin •d di + V cos jjr + V s in $ j 
d t d c 

• - . d j 

+ V cos *8 -8j + V sin — • 

( 4 ) 

The unit vectors (i, j) have an angular velocity, 


to = - >C c k, 


(5) 
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then 


= oj x i = - X r k x i = - X r j 
d t c ^ 


( 6 ) 


d J - " 


dT = w x J = 


- X c k x j = X c i. 

Substituting equations (6) and ( 7 ) into equation ( 4 ) 


( 7 ) 


r. 



- [ 

1 V cos $ - V s in -3 
L 

0 + V sin i 

3 X 

c 

i + l 


+ V cos -6 ■fl - V cos -3 X 


( 8 ) 


where equation (8) can be written as 



r 

1 


A = 

X + V sin -s X 

c 

I + 

Y - V cos i3 X 

c 


( 9 ) 


since the following terms in equation (8) are 


V cos $ - V s in $ *3 — ~ — [V cos -0 ] — -r— (X) — X. 

d t at 


( 10 ) 


V sin $ + V cos '3 -Q = — [V sin 0] = ^(Y) = Y. 


(ID 


The external force acting on the vehicle is 


i- 


f = | (F-+F cos (3-D) cos cp - N sin cp - F c sin p. sin 9 - mg cos (X c - \|f) ji 


r . 

-f j (F^+f cos(3-D)sin 9 + N cos 9 + F c sin (3 cos 9 - mg sin (X c - ''y)jj< 
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( 12 ) 



Substituting equations (9) and (12) in Newton's equation 


f = m A 


( 13 ) 


we get, 


/• 


\ 


mj X + V s in 0 | = (F^+F^ cosf3-D) cos * y - N sin 

- mg cos ( /. c - ;) 


- F c s in i sin 


(14) 


(" 

n V 


O =(F f 

- mg sin (X c - /) * 


Y - V cos j X c ) = (F^+F^ cosp-D) sin y + N cos y 4- F c sin cos ' 

05) 


At this point equations (14) and (15) can be linearized by 
using the small angle approximations, 


sin cp ^ cp sin p ^ p sin (j -f 7 ) ^ j + 7 

cos cp ^ 1 cos 6 — 1 cos ( j + 7 ) - 1 - j 7 v 1 

V- o, 

and neglecting the products of small angles* Using the above approxi- 
mations in equations (14) and (15) and solving for the highest deriva- 
tives, we obtain 


X 




- g cos X 

c 


(16) 


F^ + F 
f c 


m 


D N F 

— CD + - + — (3 + VX - g sin A . 

' m ml ^ ^ 


m 


(17) 


The degree of freedom along X is eliminated by allowing the origin of 
the coordinate system to move with the velocity Vj and acceleration 
X. The influence of this degree of freedom on the stability behavior 
of the vehicle is considered negligible. 
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The aerodynamic forces are obtained from 


D = 


C D qS 


DC 


N = 


N 


Da 


qs a = N'a 


1 ..2 1 t7 2 
q = 2 pV R ~ 2 pV * 


( 18 ) 

(19) 

( 20 ) 


as 


The thrust component normal to the vehicle x axis is defined 


r -E c e * § p - r’ P- 


( 21 ) 


Substituting equations (19) and (21) into equation (17) and denoting 
the total thrust as F = F f + F c , the resulting linearized equation of 
motion becomes 



e + v* c 


g sin X c . 


( 22 ) 


In equation (22) above the last two terms account for the 
difference between the centrifugal and gravita t iona 1 acceleration. 

Most trajectories are shaped such that the vehicle flies a "gravity 
turn ,,w trajectory. For this condition, the pitch rate X c is commanded 
such that 


g sin X 

X c = y (gravity tilt) 


(23) 


and the last two terms of equation (22) are zero. In some cases 
"gravity turn" trajectories cannot be flown due to various trajectory 
shaping constraints; in this case, the two terms should be included. 

An additional equation is required to describe the rotational 
motion of the vehicle. This equation is obtained by summing moments 
about the center of gravity as follows 


) M = Icp + d ! m + (x - x )N + (x 
/j eg Y T Cg cp eg 


XjP F c sin p = 0, 


( 24 ) 


*Also called "zero lift" and "zero angle of attack" trajectories. 



where the aerodynamic damping term d T is included. For most vehicles 
without large lifting aerodynamic surfaces, the term d 1 qp is negligible. 
Equation (24) may be written as follows using equations (19) and (21) 


.. , d' . , 

cp + — cp + 


x 


-££L 


- X 


I N'a + ■ ^ — r' p = o 


(25) 


and describes the rotational motion of the vehicle about its c.g. 

In some cases, an accelerometer is used in a control system 
to sense the accelerations normal to the vehicle longitudinal axis. 

If the accelerometer is located at the vehicle c.g., the acceleration 
sensed by the accelerometer can be obtained by summing forces along 
y, or 


r 1 

) F = m y - N - F sin 6 = 0 

u y eg c 


(26) 


and 


eg 


N' . R' 

— a + — 

m m 


(27) 


If the accelerometer is located at an arbitrary station x a 
along the vehicle, the rotational component of acceleration due to cp 
must be included or 


y 




(28) 


To control the attitude using thrust deflection, a control 
mechanism is assumed with the following form 


(3 = a Q <p + a + b Q a + g 2 y. (29) 

Equation (29) above embodies three of the more conventional 
types of control depending on the sensors used in the control system. 
The type of control is commonly called ’’gyro control” when both the b Q 
and g ; , terms are zero, ’’angle of attack control” when the g^ term is 
zero, and ”a ccelerometer control” when the b G term is zero. Generally 
both the b Q and g : terms are not used simultaneously in a control 

sys tern. 



Substituting for y in equation (28) and letting a 2 = - g2( x cg " x a ) , 
the control equation can be written in the alternate form 


P = a Q cp + a x q) + a 2 cp + b Q a + g-Ycg, 


(30) 


where 


a Q = attitude gain 

a -l = attitude rate gain 

a 2 = attitude acceleration gain 

b Q = angle of attack gain 

g 2 = lateral acceleration gain. 


The control equation shows the engine command (3 composed of signals from the 
attitude cp as measured by a space-fixed gyro, ep measured by a rate gyro, and 
either OC measured by an angle of attack meter or the local lateral accelera- 
tion measured by an accelerometer. 

In summary, the following equations are used to describe the 
pitch plane equations of motion. 


Pitch Plane Equations 


cp + 


(x - x )N' (x - x, )R' 

eg cp . eg n 


a + 




P = o 


(31) 


,'F - D\ , ( N* , 

Y = 1 ST) * + {—) a + {~ 


+ VX C - g sin X c 


(32) 


y 


eg 



(33) 


cp - 03 - 7) = a - cn. 


(34) 
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= tan 


-l 


V cos X 
W - c 

V - V sin X 
w c 


(35) 


p = a Q cp + a i9 + b 0 a + g 2 y 


(36) 


y = y cg ■ (x cg - x a> 9 


* ~ v - 7- 


(37) 

(38) 


The yaw plane equations of motion are obtained from the pitch 
plane equations by equating X=X c =y=0, y = z and Y = Z. 


Yaw Plane Equations 


(x. 


9 


_C£ 


- x )N* 
££ 


(x 


a + 


Cg_ 


V R ’ 


P = 0 



z 

eg 



p 


(39) 


(40) 


(41) 


cp - 8 = a - a 
Y w 


a 


w 


V 

w 

V 


P = a o9 + a i9 + b o a 82 'z 


z - (x 
eg eg 


- x a } <P 



(42) 

(43) 

(44) 

(45) 

(46) 
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When making fixed time point studies, the coefficients of a, j3 
and cp in the above equations are assumed constant and defined as follows: 


Ci 



x ) N' 
C P 



- V »' 

I 


Ki 



Ko = 


m 



In analyzing the control dynamics of the vehicle, it is important 
to relate the coefficients of the differential equation expressing rota- 
tional motion to the frequency and damping of the rigid body. The yaw 
plane equations given above can be reduced to the following two relations 
by eliminating a and (3 in equations (39) and (40). 


1 - g; 


(x - x ) Co + K 3 
Cg a *- 


cp + a^C^cp + 


Ci + C 2 (a Q + b 0 ) + 


+ 82 ( C 2 K 2 - C^) 


- V 


9 = 


w 


V 


Ci ~h Cobo "t g 2 ( C^Ko — C1K3) 


(47) 


(1 — goKjj) Z — a^K^cp d - a iK^cp “t 


+ Ko + (a + b Q - g 2 Ki) Kj 


9 


- (K 2 + % b 0 ) 



(48) 


Equation (47) describes the rotational behavior of the missile 
about its c.g. and is of the form 


9 + 2(; io n 9 + co^9 = F (t) . 


(49) 
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by 


Comparing equations (47) and (49), the natural frequency is given 


2 _ C -| ~t~ C P (a Q + bg) + gp(CpKp — CjK^) 

C °n 1 - go[(x - x ) C-> +* K 3 ] * 

^ eg a ^ 


and the ratio of damping to critical damping is 


i = 



(i 


a i Cp 

§2[( x cg - O C 2 + K 3 )] 


(51) 


Equations (50) and (51 r ) establish the desired relations between 
the frequency, damping and gains of the control system. 

To obtain an expression of the path reaction to the rotary 
motion the term 


Z - V 

w 

V 

from equation (47) is substituted into equation (48). The resulting 
expression is obtained: 


Z = B 0 cp + B x 4> + B^cp, 


(52) 


where 


^_a Q (CpKo - CjK^) + K 1 [C 1 + C P b n + g 2 (C P K 2 - C^)] 
c i + C 2 b 0 + g 2 (C 2 K 2 - C ± K 3 ) 


g = - " - a -J (C 2 K P - C 2 K 3 ) 

1 C-l + C 2 b 0 4 - g 2 (C 2 K 2 - C3K3) 


Ko ± b Q K P + ap(CoK P - C^) 
Ci + C 2 b Q + §2 (^2^2 " C3K3) 
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B. Drift Minimum Condition [41 


It is desirable that the gains in the control system be selected 
in such a manner that the lateral drift acceleration Z equals zero. This 
condition is called the "Drift Minimum Principle" and is determined for 
the steady state case (i.e., cp ~ q) = Z = 0). Referring to equation (52), 
then for this condition the equation reduces to 


or 


B 0 cp = 0 (53) 

- a 0 (C 2 K 2 - C x K 3 ) + KitCi + C 2 b 0 + g 2 (C 2 K 2 - C^)] = 0. (54) 
C kJ 2 ,3 ) a o = § + b o + §2 (k 2 - ^ K 3 ). (55) 


Equation (55) above establishes a relationship between the 
gains a Q , b Q and g 2 for the drift minimum condition. The drift minimum 
condition can be maintained by using angle of attack control for which 
g 2 = 0 or accelerometer control for which b Q = 0. Equations (50) and 
(51) then provide two additional equations for obtaining the gains for 
angle of attack control or accelerometer control where the frequency con 
and damping £ must be specified for the rigid body. The frequency and 
damping are selected based on the desired transient response. 

The drift minimum control case defines a control mode such that 
the sum of the force components normal to the nominal flight plane is 
zero for the steady state condition only. A claim for zero drift can- 
not be made since in transient motion cp, Cp and Z have finite values. 

During the quasi-steady-state condition the first and second 
derivatives of the variables approach zero and the attitude, angle of 
attack, and engine deflection approach limiting values. The q.s.s. 
attitude can be found from equation (47) for cp — cp = 0. 

m = f t I V A £l ± Cpbp + g g (c ? K g - C&) f56) 

^qss \ V ) C x + C 2 (a 0 + b Q ) + R; (C ; K- - C X K 3 ) ' 
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The q.s.s. values of O' and 0 are found from equations (39), (40) 

and (44) and are given below. 


a 

qss 


p 

qss 

The steady state drift in response to a wind speed V w , derived 
in Section VIII, is 


- V - 
w 


V 


J C ± + C;,(a + b ) + g : (CJC; - CjK') * 


(58) 




w 


V v 


- a. 


(C 1 + C (a + b_) + g (C :< - C.K,) 


(57) 


V 

-f - C x K 3 ) a x 

£ = 

SS C x + C.(a o + b Q ) + (C l ,K-, - ) (go + 

V 

w 

(C K_ - CbK )g + Cl + C (a Q + b Q ) ' 

1 + V L (C K; - C x K 3 ) a x j 

C. Load Minimum Control 

From equations (57) and (58) it can be seen that by reducing the 
gain a Q to zero the values of Q^ ss and P^ ss also go to zero. Since the 
forces normal to the missile axis are functions of these angles (N = N f Og 
R = F 1 (3) by reducing the gain to zero, the loads are correspondingly 
reduced. Under this condition the missile longitudinal axis is aligned 
along the relative velocity vector. Since there is no control on y, it 
can assume large values, and the missile will drift from the nominal flight 
plane. Application of load-minimum control for short durations is accept- 
able provided a deviation from the nominal flight plane is permissible. 

D. Gain Selection 

It is necessary that the rigid body natural frequency and 
damping £ be specified to compute the values of the gains. Generally, 
the selection of the rigid body control frequency depends on the fre- 
quencies of other degrees of freedom existing in the system (i.e., slosh 
modes, bending modes, etc.). In order to prevent strong coupling effects, 
rigid body control frequency should be selected below the lowest frequency 
existing in the system. A fast well damped transient response can be 
obtained by selecting a damping ratio £ between 0.4 and 0.8. 


(58a) 

(58b) 
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The gains for the control system are found from the expressions 
listed below for angle of attack control, gyro control, and accelerometer 
control. The corresponding q.s.s. angles are also listed. The gain 
relations are obtained from equations (50), (51) and (55) and the q.s.s. 
angles are obtained from equations (56) through (58). 


Angle- of -At tack Control 


a 


o 



C 2 (l + H) ; 


where 


_ £l 




H = 


(59) 


bn = H a f 


Co 


2 i u 


a i — 


(60) 

(61) 


Quasi-steady-state conditions: 


2 - V 


<P, 


w I C j ~h Cpb ^ 


qss V C x + Co(a 0 + b 0 )J 


(62) 


Z - V 


a 


w 


qss V 


— ap Cp 


Cjl + Cn(a 0 + b 0 )J 


(63) 


Z - V 


w 


qss V 


a o Ci 


Ci + C 2 (a Q 4- b Q ) 


(64) 


V 


w 


J ss 


1 + V 


fCi + Go(a o + b o ) 
(CK^ - C x K 3 ) a 2 


(64a) 


32 



Gyro Control 


a 


o 



(65) 


a i 


2 1 oj 

b n 

C :2 


( 66 ) 


Quas i-s teady-s tate angles: 


9 


qss 



V 


Cl 

Cl + C 2 a 0 J ’ 


(67) 


a 

qss 



~ a o C 

_Ci + C 2 a 0 J l 


( 68 ) 


3 


qss 



r a 


o_ 

Ci + 


Ci 

C^J’ 


(69) 


Z = V . (69a) 

ss w 
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Accelerometer Control 


GO 


- (Ci + 


S2 = 


C 2 (K X + K 2 ) - C 2 K3 + 


CgKi 

CgK g - CtK-^ 


w 2 [ (x 
n eg 


V 


C 2 + K 3 ] 


a o Kl§2 + C 2 K 2 - C3.K3 

31 = {* ' g2t(X cg " X a )C 2 + K 3 ] } • 

Quasi-s teady-sta te angles: 


9 


qss 


a 


qss 


P 


qss 


Z 

ss 


Z - V 
w 

C 1 + g P (C P Ko - C n K,) 

V 

Ci + C 2 a o + g; (C_K-; - C 1 K 3 )__ 

Z - V 


W 

- ao C P 

V 

Ct + C L) a 0 4 - go(C 2 K 4 . - C^K^) 

> 

1 

• N 

r _ 

W I 

a o 

V 

^ C j. CoS 0 § p ( Cp Kp “ C 7 K3 ) 


V 


w 


(CJC, - C 1 K 3 )g. + ( Cj_ + C :j a o ) , 
(C-Ko - C j K 3 ) a x j 



(70) 


(71) 


(72) 


(73) 


(74) 


(75) 


(75a) 
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V. FLEXIBLE BODY EQUATIONS 


A. Bending Modes and Frequencies 
1 • Homogeneous Body 
Free Vibrations of Vehicle 


The differential equation of motion of a slender beam in which 
the cross-sectional dimensions are small compared with the length and 
for which the rotary inertia and transverse shear deformations are 
neglected may be written 


m(x) 


d 2 u <j 2 

dt 5 2 k 5 


flil(x) 

I 


0" u 


w(x, t). 


(76) 


In order to solve equation (76) in closed form for the bending 
mode shapes and frequencies, consider the vehicle as a uniform beam with 
a constant mass distribution m(x) and bending stiffness distribution EI(x). 

To determine the natural mode shapes and frequencies of the 
vehicle, the forcing function w(x, t) = 0. Equation (76) can be solved 
by separation of variables such that the solution assumes the form 


u(x, t) = q(x) T’j ( t ) • 


(77) 


Substitution of (77) into (76) yields 


II 




[El (x) q(x)”]” ; 
m q(x) 


(78) 


where a prime denotes derivatives with respect to x and a dot denotes a 
derivative with respect to time. Since x and t are independent variables, 
the ratios in (78) are equated to a separation constant co 2 and we get two 
ordinary differential equations. 


t] + w 2 ij = 0 


(79) 


[El (x) q(x)"] - m(x) co 2 q(x) = 0. 


(80) 
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For a constant stiffness distribution EI(x) and 
distribution m(x) equation (80) reduces to 

constant mass 

d 4 q oj 2 n 

(80a) 

where a 2 = El/m. 


The solutions to equations (79) and (80a) are. 

respectively. 

r] = A sin wt + B cos cot 

(81) 

q(x) = C sinh \Tco7 a x + D cosh \' oj/ a ' x 

(82) 

+ E sin \j co/a' x + F cos \ ,r cj/ a 1 x. 


The constants A and B in (81) must be found from initial conditions in 
displacement and velocity at t = 0, and the constants C, D, E, and F are 
obtained by specifying four boundary conditions on the ends of the beam. 

The four boundary conditions for a vehicle unrestrained at 
each end are obtained from the condition that the shear forces and bend- 
ing moments at the nose and hinge point are equal to zero. The expres- 
sions for the bending moment and shear at any station x along the vehicl< 
axis are, respectively, 

d 2 q 

m = ex -r-t 

dx 

(83) 

and 


V , f. = EI £§ 
dx dx J 

(84) 

applying equations (83) and (84) for the conditions M = V = 0 at the 
engine hinge point, x = 0, and nose, x = j£, the four boundary conditions 
become 

0 (0) - o 0 «> ■ 0 

(85) 

( °) - o i# «> - °- 

(86) 
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Applying the boundary conditions (85) and (86) to equation (82), the 
four boundary equations are 


D - F = 0 (87) 

C - E = 0 (88) 

C cosh vWa • Jo + D sinh 4 co/ a 1 * T - E cos n/ co/ a ' • i + F sin n/ co/a ' • ^ =0 

(89) 

C sinh 4 co/a 1 • ^ + D cosh V co/a 1 • ^ - E sin >/co/ a' • Jl - F cos n/ u)/a n * £ =0. 

(90) 

The constants E and F in the above four equations can be 
eliminated to yield the two equations, 


(sinh 4 co/a ' • £ - sin nTco/ a ' • £) C + (cosh sfco/a > * £ - cos nTco/ a 1 • £) D = 0 

(91) 

(cosh \Tco/ a ' • j? - cos \Tco/ a ' • jj) C + (sinh \/co/a ' • £ + sin ^ co/a ’ • j}) D - 0* 

(92) 

A nontrivial solution to this set of equations in C and D is 
obtained by setting the determinant of their coefficients equal to zero* 
Expansion and subsequent reduction of the determinant yields the fol- 
lowing transcendental equation in terms of frequency co* 

cos 4 go/ a ’ • £ = . (93) 

cosh 4 co/a 1 • ^ 

This equation is solved graphically in Figure 13 in terms of 
the parameter \/~to7a ’ • £ 9 and it is found that an infinite set of the 
eigenvalues satisfies the equation. Also corresponding to each 

eigenvalue c^, there is a corresponding eigenfunction q n . 

From Figure 13, the frequencies (eigenvalues) for the uniform 
unrestrained vehicle are 


coq - 0 
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The frequency ojo = 0 corresponds to the rigid body mode shape. 
This can be shown by substituting co 0 = 0 in equation (80a), which will 
reduce it to 


Integrating (94), 


q = C 4 x 3 + C 3 x 2 + C 2 x + C z . (95) 


Applying boundary conditions (85) and (86) we find the constants 
C 3 = Ca = 0, and 


q = c 2 x + c l5 


(96) 


which is the rigid body mode shape. Thus, the rigid body mode shape is 
described by a rigid body translation and rotation. Generally, the mode 
shapes are defined relative to the rigid body axis. In this case, equa- 
tion (96) is unnecessary since q = 0. 

Equations (87), (88), and (91) can be solved to find the con- 

stants C, E, and F in terms of D. When this result is substituted into 
equation (82), the equation of the flexible mode shapes is 


q n oo = d 


tos n/oo /a i - cosh \' oj / a 

n n \ 


(s inh \ oj /a x + s in oj /a x) 

— # r — - — » n n 

inh \ oj / a £ - s m v oj /a y/ 


+ (coshs/oo /a" 1 x + cos \j oj /a* 1 x) j. 
n n 


(97) 


The mode shapes (eigenfunctions) q^, q^, . .., q n can be com- 

puted from equation (97) by substituting the corresponding frequencies 
oji, oj 2 , . .., ojq obtained from the transcendental equation (93). The 
mode shapes for a homogeneous vehicle with unrestrained ends (" f ree-free" 
beam) is shown below in Figure 14. The complete solution for the flexible 
motion of the vehicle with no end restraints is obtained by substituting 
equations (81) and (97) into (77) 


u(x, t) 



sin oJnt + { 


/b 


cos 


cv 


q n (V. 


(98) 
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Since equation (98) is a solution for any value n, the sum of 
the solutions is also a solution. 


u(x, t) 


y 

n=! 


A T sin (j t 4- B* cos 00 t 
n n n n 


q n W- 


( 99 ) 


Forced Motion of Flexible Vehicle 


In the previous section, the motion of the simplified vehicle 
has been found for the case of free vibrations in which no external 
forcing function w(x, t) was acting on the vehicle. For the case of 
forced motion, the mode shapes obtained previously for the free vibra- 
tion case are normalized by the deflection q n taken generally at the 
engine hinge point. In this case the mode shapes define only the 
relative displacements and will be denoted by yi(x). 


The displacement due to the load w(x, t) can be expanded in a 
series of eigenfunctions such that 


u(x, t) = \ y.(x) T) (t), 

L J 1 

n=l 


( 100 ) 


where y^(x) are the normalized natural mode shapes and T)i(t) are the 
normal displacement coordinates. Substituting equation (100) into the 
beam equation given by equation (76) yields 


00 00 

^m(x) y^x) r| i (t) + ^ [EI(x) y'Ux)]" ^(t ) = w(x, t). (101) 

n=l n=l 

The function w(x, t) denotes the forcing function acting normal to the 
rigid body center line and arises from engine forces, aerodynamic forces, 
propellant sloshing forces, etc. 


41 



Equation 101 can be reduced to a simple form by using the con- 
ditions of orthogonality between the eigen functions or mode shapes. 

The condition of orthogonality is given by 


y (x) y i (x) m(x) dx = M S 

J J -*-j 



{* 1 = J 
U i # j 


( 102 ) 


(103) 


where the eigenfunctions y^(x) and yj(x) are orthogonal to each other 
with respect to the weighting function m(x). 

Multiplying equation (80) and (101) through by y. and integrat- 
ing both sides with respect to x gives the following two equations. 


£ 


£ 


£ 


/ T| i J y i y j m dx + / T li J (EI y 'p y j 


dx = j w(x, t) y. dx 
u J 

n=l o n=l o o 


(104) 


(El y'.') y. dx 
i J 


u? / m y. y. dx. 
i J 1 3 


(105) 


Replacing the second term of equation (104) by equation (105) gives 


00 



n=l 



y i 


y j 


m dx + 


00 



n=l 



y i 


y j 


m dx 



t) y. dx, 

(106) 


which reduces to the 


fol lowing , using 


equation (102) 


r \. + oj? 

J J 



(107) 
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The quantities on the right hand side of (107) are 


£ 

M. = f my? dx 


r 

Qj = J w(x, t) y dx. 


(generalized mass) 


(generalized force) 


( 108 ) 


(109) 


Equation (107) gives the response of the jth mode to the forcing 
function Q j . If w(x, t) is independent of the motion of the beam, the 
modes are uncoupled and can be solved separately. In the actual case 
the aerodynamic forces, thrust forces and sloshing forces are coupled 
with the flexible motion of the vehicle. 


In practical applications each of the bending modes possess 
some dissipative forces which provide damping. This dissipative energy 
is small in comparison to the elastic and kinetic energies and its 
effect can be approximated by including a viscous damping term in 
equation (107). 


n. + 2t. w. tj. 

*i b i l 


2 

CJ? 71. - ~ 

1 ‘r M. 

x 


( 110 ) 


The damping ratio may vary between .0002 ^ ^ .025. 

To illustrate the response of the flexible mode to a forcing 
function assume that the forcing function is aerodynamic loading and is 
uniformly distributed over the length of the vehicle. The aerodynamic 
loading will be given by w(x, t) = w(x) sin a, where a is the angle of 
attack. If we let the angle of attack vary with the control frequency 
coc such that a = w c t, the generalized force will be 

£ 

Q_^ = sin oj Q t J w(x) y i (x) dx. (Ill) 

o 
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Then, 


w(x) y^( x ) 


T ] . + 2t . (j. Tj . + CO? T|. 
‘l 1 'l 1 'l 


M. 


sin f,j t. 
c 


( 112 ) 


The transient solution of equation (112) is obtained by setting 
Q^/M^ = 0 and is 


-too. t 

* ^ 1 . , V 

i . = A . e sin (w.t - f . ) 

‘t l d ' l 


(113) 


where oo^ is the damped natural frequency, 


= <o. \'l - . 

1 


(114) 


The particular or steady-state solution of equation (112) is 


j. = B. sin (oj t - ), 


*i l 


(115) 


where the amplitude is 


w( x ) Y t ( x ) dx 


B. = 

l 


M. oj." 

l l 


CO t= -\2. / co v 

i - -£5 ) + (2^ -^) 

u i / V <V 


;V ; 


(116). 


and the phase lag is 


2£ . — 

b i 00. 


;i = tan 


-1 


(117) 


1 - 


co. 

1 
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The complete solution of equation (112) consists of the sum of the 
transient and steady-state solutions given by equations (113) and (115), 
respectively . 


T'i . = A.e sin (co,t - \|r.) + B. sin (co t - 0 ) . 

*i i d Y i l c ^ 


(118) 


Two initial conditions are required to determine the constants 
A-l and in equation (118), and can be found by specifying that the 
initial displacement and velocity are zero; u(x, 0) = u(x, 0) = 0 at 
t = 0. The corresponding normal displacement and velocity T](0) and T)(0) 
must also be zero as can be verified by substituting these initial condi- 
tions in equation (100). The constants in equation (113) subject to the 
initial conditions q(0) = rj(0) = 0 become 



ij f. = tan* 


2 £. (1 - ^. 2 ) 1/2 




GO 


+ 1 - 




( 120 ) 


The resulting motion of the elastic displacement relative to the longi- 
tudinal rigid body axis due to the distributed aerodynamic loading is 
then 


u(x, t) 


Yy^x) 

n=l 



sin ( co , t — >]/■.) + B. sin 
d l 


(co t 
c 


0 ) 


. ( 121 ) 
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The transient term in equation (121) becomes negligibly small 
as time becomes large due to the exponential term e~£ wt . The second 
term remains sinusoidal with its amplitude and phase angle dependent 
on and co c / . The familiar resonance condition occurs as the con- 
trol frequency gj c approaches one of the bending mode frequencies oo^. 

2 . Nonhomogeneous Body 

Space vehicle^ conf igura tions are sufficiently complex in their 
structural makeup that a lumped parameter idealization of the system 
is often used in computing the mode shapes of the vehicle. The basic 
data required for the calculation of the natural mode consist of 
(1) the distributed mass m(x) and rotary inertia p(x) , and (2) the 
bending stiffness EI(x) and shear stiffness KG(x) distribution. 

For a typical space vehicle, a set of curves representing 
the structural and inertial properties are shown in Figure 15. Here 
the curves of ja(x) and KG(x) are not shown, but are of the same general 
form. In the discussion that follows, the rotary inertia p(x) and KG(x) 
will be neglected. It has been shown [1] that the rotary inertia can be 
neglected without appreciable error in mode shapes but the effects of 
shear flexibility KG(x) should be included. 

In beam theory, the total linear or angular deflection of any 
point can be represented as the sum of the deflections at that point pro- 
duced by individual forces and moments. This is a statement of the 
principle of superposition for linear systems and is expressed mathe- 
matically as 


q i 


n 




( 122 ) 


where Fj is an arbitrary force or moment called a "generalized force" 
and q-L represents the linear or angular displacement and is called a 
"generalized coordinate." The constants C^j are called the flexibility 
influence coefficients. 

If we let Fj represent forces and q^ represent linear displace- 
ment, then C^j represents the deflection at station i due to a unit load 
at station j. Using D'Alembert's principle, let the forces Fj represent 
the appropriate inertia forces; then, 


F . = - m . q . . 

J J J 


( 123 ) 
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If we consider harmonic oscillations of the system, the particu- 
lar solutions of the system are of the form 


q . = q .e 


iwt 


q . = lcj q .e 
J J 


iwt 


• * p i wt 

q . = - gj q . e 
j J 


or 


(124) 

(125) 

(126) 


F. = m. oj 2 q.e^ wt . (127) 

J J J 

Now substituting equations (127) and (124) into (122) and dividing 
through by u) 2 e lcjt , we get 


% 

to c 


C. . m . q . . 
tj J 


j=l 


When written in matrix form, equations (128) become 


91 

92 

1 

C 11 m i 

C 21 m i 

c i2 m 2 • * • c m m n 
^22 ^2 * • • m n 


9i 

9 a: 

* 

• 

• 

q n 

(a) 2 

On m i 

Cq o m 2 • • • 


_ q n_ 


(128) 


(129) 


Since the relative magnitudes of the displacements q 1; q 2 ... q n 
are sufficient to define the shape of the mode of vibration, it is con- 
venient to normalize equation (129) by the deflection q n taken generally 
at the engine hinge point. The resulting quantities will be referred to 
as components of normalized eigenvector, defined by the expressions 


yi 





(130) 
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After normalizing then with respect to q^, equations (129) become 


yi“ 




Yi ~ 

y^ 





• 

1 

“2 = 

C. . m. 


• 

• 

0J 

LJ J 


• 

y n 






( 131 ) 


Written in this matrix form, equation (131) represents a set of n 
simultaneous equations of motion and can be solved by matrix iteration 
for the natural frequency oj of the first mode of vibration and the cor- 
responding normalized amplitudes y l5 y^, . .., Yn = which define the 
shape of the elastic deflection curve* 

The square matrix on the right side of equation (131) is the 
dynamic matrix and is seen to consist of the product of the elastic 
matrix and inertia matrix, or 


C. . 


m . 
J 



(132) 


To compute the shapes and frequencies of higher modes than 
the first mode, it is necessary to suppress the dominant tendency of 
the equations of motion to converge to the shape and frequency of the 
first mode during the iteration process. This can be effectively 
accomplished by introducing the orthogonality relationship between 
principal modes of vibration, i.e., 


n 

V m y r y 8 = 0, (133) 

/ | J J J 

j=l 


where y j r and yj S are the normalized amplitude of the r and s mode, 
respectively, r ^ s # The orthogonality condition, equation (133), can 
be written in a slightly modified form by multiplying it by the ratio 
Ci n /y n r shown below to el imina te, say, the nth degree of freedom and 
yields the n equations , 


n 

V C, m. 
ln J 


V y j =0 - 


(134) 


n 


49 



The above expression can be written as follows and represents the 
modified form of the orthogonality condition. 


L L 

c in m i ^7 Yi + C in m 2 y 2 S + ... + C ln m n y n S = 0 


r r 

Vn S Vo S 

C 2 n ^ y x + m 2 y 2 + . 

y y n 

. ■'n J r 


• + C;n m n Yn * 0 


(135) 


r r 

Yi s y 0 s 

c nn m i Yi + c nn Y 2 + 

y n r y n r 


+ C nn %y n = 0 


Writing the expanded form of equation (131) to define the shape and 
frequency of the s flexural mode, we have 


S s s s 

Yi ~ co 2 (C i:l m-j_ y 1 + C 12 m 2 y 2 + • • • + Cm m n y n ) 


s s s s 

y 2 = to 2 (C 21 m 1 y x + C 22 m 2 y 2 + ... + C 2n m n y n ) (136) 


y n S = w 2 (Cni m x y x S + Cn 2 m 2 y 2 S + . . . + C nn m n y n S ). 


By subtracting the first equation of (135) from the first 
equation of (136) the second equation of (135) from the second equa- 
tion of (136) etc., y n s can be eliminated from all but the n ^ equation 
and the orthogonality relationship introduced. Performing these operations, 
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yi l c ii" c m r ) m iYi + ( C 12 -C ln r ] m 2^2 + + l c ik” c m ”) 


y 2 s =w 2 (c 21 -C 2n m iyi S + (c 22 -C 2n l2 ~) m 2 y 2 S + ... + (c^-C^ —J m k y k s 

LV y n r/ v y n r/ v y n ry 


y n s -w 2 '^C ni -C nn m iyi + ^C n2 -C nn m 2y2 ■*■•••■*" ^Cnk -< 3nn ^ 


where k = n - 1 . 


(137) 


Excluding the equation of (137), there are (n-1) equations 

in the (n-1) parameters y x s , y 2 s , •••* y n -i s - A solution to this set 
of equations is also a solution to the original n equations. 

Putting equation (137) in matrix form yields 


c ii" c in m i (c 12 -C in ^~) m 2 

y n r ^ ^ yn T// 


\ s 

Cik-Cm ] n» k yi 

y n r/ 


C 2 i-C 2 n j m i ( ^ 22 "^, 


c 2k -c 2n ) ^ ys 

y r/ 

J n 


c ki- c kn LX ~) m i 

y n 7 


c k2‘Ckn ^-j m 2 


J'k \ s 

c kk‘Ckn “) "V y k 

y n 


(138) 


where k = n - 1 and = 1 , 
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The above expression embodying the orthogonality relation can 
now be solved by matrix iteration for the natural frequency and corres- 
ponding amplitude y-^, y 2 s , . y n-1 s f° r t ^ ie s fl exu ral mode. Once 
the iteration is carried to convergence, yi s , y 2 s , y n _i s are known 

and the value of y n s can be solved for using the orthogonality relation- 
ship, equation (133), and becomes 



r s r s 

^1 yi yi + y2 y 2 + ■ • • ^ y 


k y k\ 




m y 
n 


(139) 


B. Derivation of Bending Equation [9] 

The equations necessary to describe the bending dynamics of a 
flexible vehicle are obtained by applying Lagrange 1 s equation of motion 
to the expressions for the energy and forces acting on the system. The 
energy consists of kinetic energy, potential energy, and dissipative 
energy. The energy expressions are written for three separate categories 
of the vehicle: the empty airframe, the swivel engines, and the liquid 

propellant. 


The expressions for the kinetic energy T, potential energy V, 
and dissipative energy D are given below. 

Airframe Energy 



i 


(140) 


V A = 2 


CO. 

1 


m . 
r 


1 

+ — m 
2 n 


g(x 


eg 


X af ) 


<P 


(141) 


A 



2 co 


m . 




(142) 
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where 


M* = mass per unit length of empty airframe 

I 1 = mass moment of inertia of empty airframe per unit length about 
the c.g. of elemental segment 

T| ^ = amplitude of i^ 1 bending mode 

= normalized bending deflection of i t ^ 1 mode 
= slope of bending mode 
<vu = natural frequency of bending mode 


m. 

1 


/ m 1 Y. 2 dx + / 

J 1 J 


I ! Y \ 2 dx 
o i 


£ = structural damping ratio of mode 

m f = mass per unit of length of airframe with propellant 


g 


x 


af 


F - D 
m 

location of empty airframe c.g. 


The first term of equation (140) represents the translational 
kinetic energy due to rigid body translation and rotation and also due 
to the flexible body motion. The second term represents the rotational 
kinetic energy due to rigid and flexible body motions. The first term of 
equation (141) represents the potential energy due to elastic deforma- 
tions written in terms of the kinetic energy. The second term represents 
the potential energy due to a rotational displacement in the acceleration 
field g. Equation (142) represents the energy loss of the airframe due 
to structural damping. 
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Swivel Engine Energy 


T E = 2 m E ■ (x ~~ - x > 


+ | i E ($ + ft + )’ -'I, ^(K h ) 


U 

i 


\ 1 


eg - V * + /, *‘i W " *E (* + P + £_ t \ 

i i 

(143) 


V" 


’ 2 I E “E (P ‘ + 2 “e g(x cg - V 

/ '“I ^ 

+ 8 £ E G> + 13 + ) 


(144) 


r 


“e 8 / '‘l W ^ + ) a Y -(V; 

i — -J ^ / 1 


» E * 1 h C E ( P - - Ir 5 c (ft - h'* 


E 3 E E 


(145) 


where 


I = mass moment of inertia of swivel engine about hinge point 
E 


I = mass moment of inertia of swivel engine about its c.g, 
E 


£ = length from engine hinge point to engine c.g. 

E 


“e 


= natural frequency of swivel engine 


(V = damping ratio of engine 


m E = en S ine mass - 
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Propellant Slosh Energy 


T 


s 





x o) 9 + 


- c 

Y t<*o>} 


2 



1 



- x s ) 9 


( 146 ) 




where 

m Q = fixed mass for fuel slosh anology 

I 0 = mass moment of inertia of fixed mass of propellant about its 

c.g. 


m s = slosh mass for fuel slosh anology 
x 0 = location of fixed fuel mass c.g. 

Z s = displacement of slosh mass lateral to vehicle centerline. 
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x s = location of sloshing mass c.g. 

ojg = natural frequency of oscillating propellant 

t = damping ratio of propellant. 

In order to use Lagrange’s equation, the following derivatives 
are necessary. 



1 


(149) 




Y. 

i 


dx 


9 


m 


A af 


Y. 

i 


dx + 


V 

L 


r *i 


- $ 


I' Y’ dx 
o 1 



I' Y'. 2 dx 
o 1 


m' Y? dx 
A l 


(150) 



(151) 




i 


(152) 
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0 


( 156 ) 



where Sg = £%. 



( 157 ) 


+ ) m s ^Z-^ + Y.(x s ) + Z s Y.(x s )j 


. c)T 
d s 


ie — L = Y i (x o ) 2, ra ° t z ■ + 23 Y i (Xo> ; 




- Y^(x 0 ) } I 0 ^cp +^1). Y'(x, 




+ Y.(x g ) ) m s jZ - £ g cp + ) I], Y.(x g ) + Z s f 


bV r-i 

= Y'. (x ) x / m g Z 

^T) i i s s 6 s 


( 158 ) 


( 159 ) 


dD 

= 0 . ( 160 ) 

Using the above derivatives and the Lagrange equation, the 
bending equation becomes 
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m' Y . dx + ) m n Y . (x ) + ) m Y . (x ) - Z + 
A 1 / . ° 1 o / . s i s / 1 


m' £ , Y. dx 
A af l 


h Y i (x h } “ / m o^oYi(x 0 ) - > I D Yl (x Q ) - 


m £ Y. (x ) 

S S 1 s 7 


X o Y f dx ^ + , Y i 


k dx - 


r y' dx - I E r. (Xh ) 


+ Y ?(x„) - ^ I„ Y^(x 0 ) + 


I", Y ?<vb- h w p 


+ )m Y.(x ) Z + ) gj . 2 m. ti. + ) 2to. m. t. fi. 
/_, s i s s s Z_i 1 11 /_, 1 1 b i 'i 


+ "E W 1 Z - V + /.‘ T i W $ + * + ),\ Y *i<V 


11 

/ 


' S E Y i ( V { £ - t cg'f + /Ji W} + r E Y i (x h ) + 
tS I* Y 'i ( V {<P + Y 'i (x h ) } -"e i Y 'i<V 1 11 i W 


”e i Y i ( V f ± ’ll Y i 



S 1 

\ 


+ g % Y s - q„ • <i6i) 


Equation (161) can be simplified by noting that the coef- 
ficients of Z and cp represent mass and inertia terms taken about the 
vehicle c # g. Since the linear and angular momentum is conserved for 
the normal modes of vibration, then 



Y.(x) 


dx + 


\ m Y. (x ) + 

L 0 i o 



Y. (x ) = 0 
l s 


l 


( 162 ) 
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and 


v 


m I *af Y i (K) dX + l E Y i (X E } " ^ m o^o Y 1 (x q ) - ^I 0 Y’ (x Q ) 


+ ) m H Y.(x ) + / I* Y. (x) dx = 0 
7 ss s' I o 


(163) 


and, therefore, the first two terms containing Z and cp drop out. 

The coefficient of T] in equation (161) can be defined as the 
generalized mass by letting 


m. = / m' Y 2 dx - / I' Y'. 2 dx - !L Y'. 2 dx + 

1 ,/ A 1 , oi Ei 


m 0 Y? (x Q ) 


Io y;. 2 (x 0 ) + V m Y? (x s ) . 


C_J 


S 1 


(164) 


Now using expressions (162), (163) and (164) and regrouping 
some of the engine terms, the bending equation can be written as follows 


m.q + 2m. to. f . n. + m. or n . + ) m Y.(x ) Z - g ) Y 1 . (x Q ) m Q Z s 

X ' 1 X b X ‘x 1 l'l / . s i s s & / . x s . s b 


- X E n (x h> P + -E W - S E W z - v + 7 . v ‘i W 


- VeW - t e w m 5 + p + / ,‘ T ii Y, <v 


”e i {w I 7 *i w + ^ n ( v [ 9 + p+ w 


- W <p + MW 


= Q. 
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V 


(165) 



The generalized forcing function Q T ^ is defined as 



t) Y^(x) dx 


(166) 


where w(x, t) is the external forces acting on the vehicle and Y-^(x) 
defines the mode shape. The two primary forces which affect bending 
are the thrust and aerodynamic normal forces. They are as follows: 


Thrust 


w T (x, t) « R'0 


Q (T) = [R' Y.(x h )] (3. 

Tl ± 1 11 


(167) 

(168) 


Aerodynamic 


a 2 


w A (x, t) 


= (x) . a q S; 


KX 




s 


a 


ox oCt 


(169) 


<0 - q S 
1 


cjj ( x)a Y i (x) dx 


= q s ) cfj (X n ) a (Ax n ) Y i (x n ). 
L i 

n 


(170) 


then 


The generalized forcing function of the bending equation is 


Q = 

T 'i T ‘i T ‘i 

[R' Y.(x h )] 0 + q s J^Ch a (x n ) a (Ax n ) Y.(x n ). 

n 


(171) 
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The bending equation can now be written as 


H. + 2f .Uj.fi. + LK 71. = 

‘l 1 'l L 'l 


( 172 ) 


where 


W v .. g w v 1 

in. / m sk ^sk + m. /_, m sk Z sk 

1 k 1 k 


Y ,(x u ) I Y (x, ) - « Y' (x ) + I Y' (x, ) Y' (x, ) 


r L V* [Y x h 

j 


e vvj + ^ vv yv 


m E g ^ r 


s e w YVj V — 4 tyv Y j (x h } + W w 


*E W Y j ( Vj T »J - ^ L Y i (X h> - ^E WJ z 


(173) 


x - x, + g ) Y.(x ) - S (x - x.) Y'. (x, ) - I Y'. (x, ) r 9 

eg h E i h E eg h l h E i h j x 


+ (y. (xj - g v Y'. (xjjtp + ^ { [ c' (x)aY.(x) dxj 


[ i h E h'J Y m. 


r { S E W - X E w} e • m 1 l s E i W • R ' Y i ( V) p - 
1 ^ J 1 ^ ^ 


Equation (172) above describes the bending dynamics of the flexible 
vehicle for each i t ^ 1 mode of vibration. Equation (173) represents the 
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forcing function due to k propellant slosh modes, engine dynamics and 
aerodynamic forces. Frequently, the r]j , T) . , Z, cp and cp terms are small 
in comparison to the other terms in equation (173) and are neghected. 

The open-loop coupling between each orthogonal bending mode 
(t] ^ to ij j ) can occur through the following forces: (1) aerodynamic 

damping, (2) engine thrust and engine inertia, and (3) engine gimbal 
bearing friction. Some of these effects have been neglected in the 
forcing function of equation (173). 


The open-loop couplings between the nonorthogonal modes in equa- 
tion (173) are 

(1) T] to mode Z g , coupled inertia 1 ly and elastically, 

(2) T) ^ to mode p, coupled inertially and through thrust forces, 

(3) to mode a, coupled through aerodynamic forces. 


The closed loop coupling of r \ ^ to mode tj. depends on the trans- 
mittance from mode to control system, i.e., the way in which the mode is 
sensed by the control system. 


VI. SLOSH EQUATIONS [9] 

The motion of liquid propellants in cylindrical tanks can be 
described by the well known spring-mass model. This method gives the 
engineer an insight into the problem and lends itself to a description of 
the liquid motion which is readily adapted for digital and analog simulation. 

This mechanical analog duplicates the hydrodynamic forces and moments 
acting on the tank which are simulated by replacing the slosh mass by an 
equivalent rigid mass for each fundamental slosh mode. Reference 2 con- 
tains the necessary graphs to determine the slosh mass, slosh mass loca- 
tion, associated spring constant and damping constant, and other slosh 
variables versus the liquid column height- to- tank-d iame ter ratio for 
tanks with circular and annular cross sections. 

The balance of the propellant exclusive of the slosh mass is assumed 
rigid and is included in the total mass and mass moment of inertia of the 
vehicle. 

The pendulum analogy used to duplicate the motion of the sloshing 
propellants is similar to the spring -mass model. In this case, the slosh 
mass is replaced with a rigid mass suspended by a massless link in length 
£ as shown below. 
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PENDULUM SLOSH MODEL 


In this case, the length £ is chosen in such a way that the ratio 
of the longitudinal acceleration g to the length of the link £ represents 
the square of the natural frequency of the liquid mode. This model 
describes the motion of the liquid propellants equally as well as the 
spring mass model, but the equations will be shown for the spring mass 
model only. 

Generally only the fundamental mode of oscillation is considered 
since the slosh masses of successively higher modes have less influence. 
Also test experience indicates that a great deal of turbulent mixing 
accompanies high frequency tank oscillations, and that damping effects 
are greater, further reducing their significance in control analysis. 

The slosh model and the flexible body configuration are shown in 
Figure 16. The equations of motion for the slosh mass using the spring- 
mass-damper system can be conveniently derived using Lagrange 1 s equa- 
tions. Only the fundamental slosh mode will be considered for each of 
ff j ,r tanks mounted along the vehicle center line. It will be further 
assumed that the motion of the slosh mass is only influenced by rigid 
body translation, rigid body rotation, and bending oscil la t ions . 

Before applying Lagrange’s equation an expression for the kinetic 
energy T, potential energy V and dissipative energy D must be written. 

^ 2 

T = \ m s |z - i s q> + £ t,. Y.(x s ) + Z g j (174) 

i 
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FIGURE 16. SPRING MASS DAMPER SLOSH MODEL 
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\ k S z s + m s 5 ^s'P 2 " m s i z s {9 + Y i (x s ) 


(175) 


4 c 2 2 . 

2 s s 


(176) 


The Lagrange equation is 


A A- , 5v . Sd n 

dt <$q. dq. bq . ^i 5 

1 L 1 


(177) 


where the generalized coordinate, q-j_ = Z s , is the one of interest. 
Taking the respective derivatives, we get 


If = m s X + i s ' £ a$ + 1 *>1 Y i (x s } 


(178) 


■hw s -\i i + K ■ W 

, £ 


(179) 


r r — 1 


= k Z - m g ^ cp + ) q . Y* (x ) 
dZ s s s 6 l 'i 1 s 


(180) 


= c Z 
dz s s 


(181) 


Q = 0. 
z 


(182) 


The slosh equation now becomes 


m JZ + Z 
si s 


- ; \ z s - m s 8 9 + > t), n(- s )r + v = 


(183) 
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Since 2£ s ca s = C s /m s and w| = k s /m s , the slosh equation becomes 


Z + 2t a) Z + w 2 Z 
s ^s s s s s 




Y^(x 




i 


(184) 


The right-hand side of equation (184) represents the forcing 
function of the sloshing mode. The propellant sloshing modes are forced 
only by inertial coupling with the vehicle rigid and elastic body modes. 
The coupling of the sloshing modes back into the vehicle is a combination 
of both inertial and elastic coupling. 


If the flight control system directly senses propellant motions, 
then it responds to these signals and reacts on the body modes; these 
in turn affect the propellant modes. If the flight control system does 
not sense propellant motion directly (the more common case), then it 
sees these motions only as the body modes (which is the case with con- 
ventional sensors). In either case, the slosh modes must have their 
loops "filtered" by the body modes. It follows then that these pro- 
pellant modes will only couple significantly with modes in their 
immediate frequency range or with modes having a large enough pass 
band to encompass the slosh frequencies. 


VII. WIND REPRESENTATION 

The following is a brief outline of the method to be used in 
applying the MSFC "synthetic wind profile" to control system studies. 

The quasi-steady state nondirectiona 1 wind condition used will not 
be exceeded more than 5 percent of the time for the worst monthly period. 
This corresponds to a 95 percent wind speed profile envelope as shown in 
Figure 17. For example, from this profile the value of the wind speed 
is 75 meters/seconds in the 10 to 14 km region, commonly referred to as 
the "jet stteam" region. The aerodynamic loading is typically most severe 
in this region since this is also in the range of altitudes where the 
maximum dynamic pressure occurs during the boost phase of flight. 

Figure 18 gives the 99 percent probabil ity-of-occurrence vertical 
wind speed change spectrum envelopes as a function of altitude (y) and 
scale-of-dis tance (Ay). These curves define the x^ind buildup rates 
(shears) which for the given scale of distance will not be exceeded 
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more than one percent of the time for the worst wind month. The shear 
data arecompiled for altitude layers to 5 km in depth. By starting with 
the wind speed given by the profile envelope for the altitude of interest 
and subtracting shear for each altitude layer, the wind buildup may be 
determined over a 5 km layer. The procedure for obtaining the wind buildup 
is shown in Table I where the selected altitude was 13 km. The 95 percent 
wind speed corresponding to this altitude is obtained from Figure 17 and 
found to be 75 m/s. The values of AY and AV W shown in Table I are then 
obtained from Figure 18. These tabular data are then plotted as shown in 
Figure 19. 

The angle of attack due to wind, 0 ^, may now be constructed by 
using the relation 


or 


V cos X 

a = 77" — 7 : 7 (pitch plane) 

w V - V sxn X 
w c 


V 

a = ~ (yaw plane) 

w V 


where V is the vehicle velocity in (m/s) and X c is the tilt angle measured 
from vertical at launch. The Qfo curve can be computed for a frozen 
time-point study (constant coefficient study) or for a time-varying 
s tudy . 

For a constant study, the V and X c at t = t [ (Vw max ^ ] could be 
used or the velocity, and tilt angle could be used corresponding to the 
trajectory time from t x = t[h (V WmaK ) - 5km] to t P = t[h (Vw^^) ] . The 
wind curve V Xv r(t) appears in Figure 20. A line can be faired into the 
curve at t x from V w (t) = 0, without affecting peak transient values. 

The wind starting point should be 2 seconds before t 1# The wind max 
plateau (P) may hold for a short time or simply be a peak. There are 
no published data which state statistically how long the wind maximum 
conditions may exist. From past observations, the altitude layer over 
which the maximum wind peak may be held at the 95 percent maximum wind 
speeds will not exceed three kilometers. After the plateau, wind shears 
^ 99 percent shears may be used to fair the curve into the quasi-steady- 
state wind speed envelope curve. 

The relationship between established gust and/or embedded jet 
characteristics and the idealized wind speed profile envelope is shown 
in Figure 21. The embedded jet is found in the 4 to 15 km region, and 
is statistically independent of the wind shears. The gust in Figure 21 
is a typical one. The gust may be superimposed on the top of the 
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quasi-steady-state wind curve in Figure 20. When the gust and 99 percent 
shears are combined in this fashion, they should be reduced by 15 percent 
to maintain a more likely probability level. This method is applied in 
Table II. 

Figures 22 and 23 illustrate the typical response for a large 
vehicle for two types of wind disturbances. In the first example, the 
curve contains the 99 percent shear buildup into the 95 percent pro- 
file. In the second case, a quasi-square wave gust is added. 

The preceding figures are used primarily for a rigid vehicle. For 
nonrigid studies, sinusoidal gusts are of major interest. For simula- 
tion of the sinusoidal gusts, wind speed changes for Ay scale distances 
< 1000 meters should not be used to construct x^ind buildup rate portion 
of profile. The gust \\rave length, the peak-to-peak amplitude, and 
number of successive gusts are given in Table III for use with the 
quasi-steady-state winds. The sinusoidal gust curve will appear as 
shown in Figure 24. 

For a time-varying x^ind response study starting at liftoff, the 
synthetic x^ind profile constructed in Figure 19 is used except that the 
curve will originate from an altitude of zero and connect x^ith the curve 
representing the statistically designed wind shear (x^ind buildup rate) 
curve at 8 kilometers altitude as shox^ in Figure 25. To obtain the 
\tfind relationship above this altitude (13 km), one may back off from 
the peak using the same wind speed change rates as for the buildup, and 
then fair into the quasi-steady-state wind speed envelope curve. Also 
one could fair the curve to the x^ind speed envelope curve with less (AV W ) 
xtfind speed change rates than for the buildup curve. The wind curves 
should be constructed for any altitude in the same fashion. 

The wind tables and graphs included in this section are shown ‘for 
illustrative purposes only. Since wind data is constantly being revised 
and updated, it is suggested that the latest wind data be obtained from 
the Aerospace Environment Office, Aero-As trodynamics Laboratory, before 
attempting to apply the methods for xvind construction discussed in this 
section . 
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TABLE I 


WIND SPEED BUILDUP FOR 8 KM TO 13KM ALTITUDE 


Ay 

Meters 

y 

Meters 

95% Wind Speed 
at 13 km 
m / sec 

AV W 

for Altitude Layer 
m/sec 

Wind Speed 
m/ sec 

0 

13000 

75 

0 

75 

100 

12900 


9 

66 

200 

12800 


14 

61 

400 

12600 


22 

53 

600 

12400 


27.4 

47.6 

800 

12200 


31.7 

43.3 

1000 

12000 


35.3 

39.7 

2000 

11000 


48 

27 

3000 

10000 


55.2 

19.8 

4000 

9000 


59.2 

15.8 

5000 

8000 


61.5 

13.5 
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TABLE II 


WIND SPEED vs ALTITUDE INCLUDING GUST 


Ay 

Meters 



y 

Meters 

95% Wind Speed 
at 13 km 
m/sec 

AV W 

reduced by 15% 
m / sec 

Wind Speed 
m/sec 

+ 

50 -> 300 

13050 -> 13300 


0 

75.0 

+ 

0 250 

13025 -> 13275 


7.65 

82.65 

+ 

25 

13025 


7.65 

82.65 


0 

13000 

75. 

0 

75. 


100 

12900 


7.65 

67.35 


200 

12800 


11.9 

63.1 


400 

12600 


18.7 

56.3 


600 

12400 


23.29 

51.71 


800 

12200 


26.945 

48.0555 


1000 

12000 


30.005 

44.995 


2000 

11000 


40.80 

34.20 


3000 

10000 


46.92 

28.08 


4000 

9000 


50.32 

24.68 


5000 

8000 


52.275 

22.725 
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TABLE III 


SINUSOIDAL GUST CHARACTERISTICS [13] 


A, Gust Wave Length 
(meters ) 

A, Peak-to-Peak 
(Amplitude (m/s) 

No. of Successive 
Cycles 

50 



3.5 

10.0 

100 

5.0 

8.0 

200 

9.0 

6 . 0 

300 

11.5 

5.0 

400 

13.5 

4.0 

500 

15.0 

3.0 

600 

16.5 

3.0 

700 

17.5 

2.5 

800 

18.5 

2.0 

1200 

21.0 

1.8 

1600 

23.0 

1.7 

2000 

24.8 

1.6 

2400 

26.0 

1.5 
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Altitude 


Scale-of-Distarjce (Distance below given altitude 


J 

s 

s 

s 

<1 

o 

o 

o 

o 

o 

o 

O' 

o 

o 

o 

r-H 

CM 

on 

<1- 



To account for the effect of wi d directio 
reversals, the 20 to 35 km altitude wi d 
speed change values have been increased 
by 25 from those given in report, 

' ASA T.i D-1274 


Wind Speed Change, AV W (m/sec) 


FIGURE 18. NINETY-NINE PERCENT PROBABILITY-OF-OCCURRENCE VERTICAL WIND 
SPEED CHANGE SPECTRUM ENVELOPES AS A FUNCTION OF ALTITUDE 
AND SCALE-OF-DISTANCE ASSOCIATED WITH THE NINETY-FIVE AND 
NINETY-NINE PERCENT WIND SPEED PROFILE ENVELOPES FOR 
ATLANTIC MISSILE RANGE [11] 
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<TLD SPEED BUILDUP FOR 8 KM TC 13 ALTITUDE 
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FIGURE 21. RELATIONSHIP BETWEEN ESTABLISHED GUSTS AND/OR EMBEDDED 
JET CHARACTERISTICS (QUASI- SQUARE WAVE SHAPE) AND THE 
IDEALIZED WIND SPEED PROFILE ENVELOPE 
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FIGURE 22. TYPICAL VEHICLE RESPONSE TO A 
WIND BUILDUP WITH MAXIMUM PLATEAU 



FIGURE 23. TYPICAL VEHICLE RESPONSE TO 
WIND BUILDUP WITH GUST 




FIGURE 24. WIND SPEED vs TIME FROM 0 KM TO 13 KM 




VIII. TRANSFER FUNCTIONS AND BLOCK DIAGRAMS 


The transfer functions contained in the block diagrams are most 
generally shown in the Laplace variable S. 

The Laplace transform of a time function f(t) is defined as 


F (s ) = L[f (t) ] = J f(t) e" St dt, 
0 + 


where S is a complex variable. For F(s) to be meaningful, the integral 
must converge and f(t) must be defined for t > 0 and equal to zero for 
t ^ 0. There are limitations on the functions f(t) that are Laplace 
transformable. To assure convergence of the Laplace integral requires 
that (1) f(t) be piecewise continuous over every finite interval and (2) 
that f(t) be of exponential order. 

Consider for example, the Laplace transform of the derivative of x(t). 

It is 


dx 

dt 


r dx -st 

■ J 51 e dt ’ 

0 + 


which gives the result 


dx 

dt 


= S X(x) - x(0 + ) , 


where x(0 + ) is the initial condition. 


Similarly, the Laplace transform of the second derivative of x(t) is 


L 

V 2 


= S 2 X(s) - S x(0 + ) - d 1 
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where 


s x(0 + ) 


and 

dx(0 + ) 

dt 

are initial conditions. The Laplace transform of other functions commonly 
used in control theory can be found in References 7 and 8. The outstand- 
ing features of the Laplace transform are its ability to convert a dif- 
ferential equation into an algebraic equation, its ability to handle 
initial conditions directly, and its wide range of application. 

A transfer function expresses the dynamic relationship between out- 
put and input. In general, the output quantity of any linear component of 
a system is related to the input by a gain factor and combinations of 
derivatives and integrals with respect to time. 

For example, consider the second order differential equation 


x + ax + bx = K(y + cy) . 


The right-hand side represents the input or forcing function of the system, 

K represents a gain factor, and x represents the output or response of the 
sys tern. 

The Laplace transform of this equation is 
S 2 X(s ) - S x(0 + ) - dX d ( ^ ' + a [SX(s ) - x(0+)] + b X(s) = K[SY(s) - y(0+) + CY(s)] 
and 

. Kjs ± C) Y(S) + s H- ^ +,■(<>+) -y(0+> 

K J S 2 + aS + b S 2 + aS + b 

The first term on the right-hand side of this equation represents 
the transfer function and the second term represents the initial condition 
operator. Frequently in control system studies, the initial conditions 
are made zero for simplicity. In this case for zero initial conditions, 
the ratio of output to input becomes 


X(S) _ K(S + C) 

Y (S ) S 2 + aS + b‘ 
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The roots of the numerator in the above expression are called zeros, and 
the roots of the denominator are referred to as poles. 

The Laplace transform of the yaw plane equations shown in Section IV A 
are summarized below assuming zero initial conditions. These equations 
are written in terms of the wind angle C^(s) and represent the rigid body 
equations using a control law of the form 


P = a D cp + + b 0 a + g 2 z. 


A(s) = (1 - g 2 K 3 ) S 3 + 


i (K 2 + b 0 K 3 ) + ai C; : S l 


Ci + C.,(a 0 + b 0 ) + (C 2 K 2 - CiK 3 ) + g 2 ) j S (184a) 

_ 


1 

+ V 


(C : ,K ; . - C 1 K 3 )(a c - K 1§2 ) - (C x + C,b 0 ) K x 


Cp(s) = - 


OC (s) 

w 

A(s) 


(C 2 K 2 - CjK 3 ) g 2 + C 2 b + C- 


(184b) 


a(s) = 


% (a) 

A(s ) 


(1 - K^go) “f - CAa 1 S C->a 0 J S 


(184c) 


(3(s) = 


a (s) 

W 

A(s) 


(b Q + K 2 g 2 ) S 2 - CjajS - C iao S 


(184d) 


Z (s ) 


a (s) 

w 


(K 2 + K 3 b Q ) S 2 + ( C 2 K 2 - C x K 3 ) a ± S 


+ (C 2 K 2 - C^K 3 ) (a Q - K x g 2 ) - (Ci + C- b D ) Ki 


(184e) 
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(K 2 + K 3 b G ) S 2 + (C 2 K 2 - C X K 3 ) a x S 


z(s) = 


a (s) r 
A(s) 


+ (C 2 K 2 - C 1 K 3 ) a Q 


S. 


(184f ) 


Equation (184a) above is called the characteristic equation of the 
system. Equations (184b) through (184f) are transfer functions of the 
respective variables with respect to the wind angle C^(s). When drift 
minimum gains (a 0 , b Q , g 2 ) are used, the constant term in equation (184a) 
and the drift equation (184e) are zero, or 


(C 2 K 2 - CJC 3 ) (a 0 - K lg2 ) - (C 2 + C 2 b 0 ) K x = 0. (184g) 


In this case, an S can be factored out of the A(s) equation and the 
numerator of the 2 (s) equation as follows: 


a(s) dm = ^K 1 - S 2 K 3 ) S 2 + 


-(K 2 + b Q K 3 ) + a X C 2 


S + 


+ (C 2 K 2 - C x K 3 ) (-1 + g 2 ) 


Cl + C 2 (a 0 + b Q ) 

(184h) 


Z (s ) = 


a (s) 
w ' 7 


A(s) 


DM 


(K 2 + K 3 b 0 ) S + (C 2 K 2 - C x K 3 ) a ± 


S, 


(184i) 


where the subscript DM denotes M drift minimum . 11 

Notice that substituting A(s)^ in equations (184b), (184c), (184d), 
(184f) and (184i) , the S factor in both the numerator and denominator cancel 
as a result of the drift minimum gain relation (184g) . Solving for the 
three roots of A(s) D ^ = 0 in equation (184h) reveals that one root is located 
at S = 0. Since this root is located at' the origin of the S plane, it is 
commonly referred to as the "drift root pole." 
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An additional Laplace transform theorem which is useful in control 
theory applications is the final value theorem . This theorem states that 
if a function f(t) and its first derivative f* (t) are Laplace transformable 
and if the poles of S F(s) lie ins ide the left-hand S plane, then 


lim S F(s) = lim f(t). 

S-» 0 t— > oo 


This theorem is applied 
the steady state drift for a 


below to the expression for Z(s) to find 
wind profile given by 


a (t) 

W 




where V w /V is a constant. The transform for this expression is 


(V /V) A 

Q W (S) = S(S + Ai • 


(184 j ) 


Substituting (184h) and (184j) into equation (184i) and applying the 
final value theorem gives 


1 im Z ( t ) = 

t— »oo 


1 im 

0 


r<v w /v)A 
S(S + A) 


S[(K 2 + K 3 bo) S + (C 2 K 2 - C 3 K 3 ) aj 


r 

JLS{(l-g 2 K 3 )S 1 2 +[i(K 2 +K 3 b o +a 1 C 2 ) ]S+C 1 +C 2 (a 0 +b 0 )+(C 2 K 2 -C 1 K 3 ) (^hgO }- 


a nd 


(VjV) (C 2 K 2 - C x K 3 ) a x 


SS C x + C 2 (a 0 + b G ) + (C 2 K 2 - C x K 3 )(2i + g 2 ) 


(184k) 
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The same result could have been obtained for Z ss by using a step 
wind Qfo(t) = (V w /V) or any other wind profile which converges to a final 
value (V w /V) as t o>. Hence, the -final steady state drift is independent 
of the wind profile Qfo(t) as long as it converges to a final constant 
value. 


To illustrate application of Laplace transform theory to obtain 
transfer functions and the subsequent construction of the block diagram, 
consider the following simplified rigid body equations with attitude (ep) 
and attitude rate (q;) ^control. The drift (2) and wind angle (o^) are 
ignored in these equations. 


cp + C-jQ: + C 2 p = 0 
cp - a = 0 

P = a o(q? - Cp c ) + a'xcp. 


In this system of equations, cp c is the commanded attitude of the vehicle 
and acts as the forcing function. The output of the system or the con- 
trolled variable is cp* 

Taking the Laplace transform for zero initial conditions, solving 
for the transfer function cp to p and using the first two equations results 
in the following rigid body transfer function 


<p(S) = - Co 

p(S) S 2 + C-L * 

This can be represented by the following block diagram and is referred 
to as open loop control . 


p(s) m 

~ c g. 

Cp(s) 


s 2 + c x 
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The third equa tion in the above sys tern of equa tions represents the 
engine deflection 3 in response to the commanded attitude cp c and feed- 
back from the attitude gyro cp and rate gyro q>. This equation is repre- 
sented by the following blocks and summing junction. 



a 

o 


c p c (s) 


Combining the above two block diagrams results in the following closed 
loop attitude control system. 


cp^ s ; 

I a ] 



1 ° 1 





i 3(a) 


■RIGID BODY 


C-.> 


cp(s; 


u 


a o + 3 iS 



POST '": A. D 

-'ATP JiAO 


Host practical control systems are much more complicated than the 
simplified illustrative example given above and often contain multiple 
feedback loops and several inputs. 

By means of block diagram reduction theorems (References 7 and 8 ), 
every block diagram can be reduced to isolate parameters in a transfer 
function (if necessary to study its effects) or can be reduced to the 
basic block diagram as shown in Figure 26. The block diagram reduction 
provides a technique to reduce the system to a form which can be studied 
in terms of general feedback theory. 
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The block diagram shows a pictorial representation of the control 
systems operation and the manner in which signal information flows 
throughout the system. A rectangular block represents a transfer 
function, lines represent paths of signal flow, a rrows indicate direc- 
tion of signal flow, and circles represent summing junctions. 

The basic block diagram of a negative feedback control system is 
shown in Figure 26. 



FIGURE 26. BASIC BLOCK DIAGRAM FOR A FEEDBACK SYSTEM 


The variables 
R(s) 
E (s ) 
B (s ) 
G(s) 
H(s) 
C(s) 


in the block diagram are defined as follows: 
= reference input to control system 
= error signal 
= primary feedback signal 
= forward transfer function 
= feedback transfer function 
= controlled variable. 
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From the block diagram shown above, the following relations can 
be written relating the variables of the system. 


E (s) = R(s) 

- B (s) 

(185) 

C(s) = G (s) 

E (s) 

(186) 

B(s) = H(s) 

C(s). 

(187) 


Solving equations (185), (186) and (187) for the control ratio 
C(s)/R(s), we obtain the closed-loop transfer function 


c O) = G ( s ) 

R(s) 1 + H(s) G (s) * 


(188) 


Solving equations (186) and (187) for B(s)/E(s), we obtain the 
open-loop transfer function 


ten - h < s > g < s >- 


(189) 


The characteristic equation of the system is the denominator of 
the closed-loop transfer function or 


1 + H (s) G (s) = 0. 


(190) 


The foregoing equations are written assuming both G(s) and H(s) are 
positive. Obviously, G(s) and H(s) may occur with other sign combina- 
tions. Negative feedback will always occur, though, when the feedback 
signal B(s) and input signal R(s) differ in sign, and positive feedback 
will occur when they both agree in sign. Positive feedback systems are 
often unstable since this signal has the correct phase for regeneration 
to take place around the feedback loop; and, if the amplitude of the 
feedback is sufficient, sustained oscillations will result. 
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IX. METHODS OF STABILITY ANALYSIS 


A. Routh 1 s Stability Criterion 

Routh's Stability Criterion is a technique for determining the 
number of roots of a polynomial equation with positive real parts. 

When applied to the characteristic equation of a control system, (e.g., 
equation (190)), it provides a simple test for absolute stability 
because a stable control system cannot have any roots of the charac- 
teristic equation with positive real parts. 

Let the following equation represent the characteristic equa- 
tion of a control system: 


B n S n + B n _ jlS 11-1 + B n . 2 S n " 2 + ... + B 0 = 0. 


(191) 


The coefficients of the characteristic equation are arranged in the 
first two rows as shown below. These coefficients are then used to 
evaluate the rest of the constants to complete the array. 


s n- 


s n-2 


s n- 


B~ 


Bn-2 ^n-4 B n _ 6 


Bn- 1 B n-3 B n _ 5 B n _ 


Ci C z 


di d- 


(191a) 


The constants C l9 Co, C 3 , etc., are evaluated as follows; 


Ci = 


Bn- 1 Bn-s ~ Bn Bn-; 
Bn- i 


(192) 
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Co = 


B n -1 B n-4 “ B n ^n -5 


Bn- 


n-i 


(193) 


Bn-i B n «6 - B n B n _ 7 


C-, - 


B n-i 


(194) 


This pattern is continued until the rest of the C' s are zero. The 
constants d l3 d 2 , d 3 , etc., are formed using the S 11 " 1 and S n “ 2 row and 
are 


c i B n _ 3 - B n _ 2. C 2 



(195) 


c i B n _ 5 - B n _ j_ C 3 



(196) 


This procedure is continued down to the S 1 and S° rows, which 
contain one term each. Once the array has been found, Routh T s criterion 
states: The number of roots of the characteristic equation with positive 

real parts is equal to the number of changes of sign of the coefficients 
in the first column. Therefore, the system is absolutely stable if the 
terms in the first column have the same sign. 

In forming the Routhian array, the following theorems may be 

used: 


1. The coefficients of any row may be multiplied or 
divided by a positive number without changing the 
signs of the first column. 

2. When the first term in a row is zero but not all 

the other terms are zero, substitute a small positive 
number 5 for the zero and proceed to evaluate the 
rest of the array; or substitute in the original 
equation S = 1/x and then solve for the roots of x 
with positive real parts. The number of roots x with 
positive real parts will be the same as the number of 
S roots with positive real parts. 
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3. When all the coefficients of one row are zero, 
the Routhian array can be completed by replac- 
ing the zero row with the coefficients obtained 
by differentiating the auxiliary equation 
formed from the preceding row. For example, 
assume all the coefficients d ±y d^, ... of 

the S n ~ 3 row of (191a) are zero; the auxiliary 
equation obtained from the S n ~ 2 row is 


c x S n ~ 2 + C 2 S n ~ 4 + C 3 S n "° + . . . = 0. 


To complete the array, the auxiliary equation 
is differentiated and becomes 


(n - 2)C X S n “ 3 + (n - 4)C 2 S n " 5 + (n - 6) C 3 S n ~ 7 + ... = 0. 


The coefficients of this equation are inserted in the 
S n ~ 3 row of equation (191a) and are 


d 1 - (n - 2)Cj 
d 2 = (n - 4) C P 


In order to illustrate an application of Routh’s criterion, 
consider the following system of equations for a rigid vehicle with 
attitude, attitude rate, and angle of attack control: 
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cp + c x a + c 2 p = 0 (197) 

Z - K lC p - - K 3 P = 0 (198) 

| = cp - a + (199) 

P = a 0 cp + a iq) + b Q a. (200) 


This system of equations can be reduced to three equations 
in three unknowns by eliminating the Z variable between equations (198) 
and (199). Then taking the Laplace transform of the system for zero 
initial conditions of the variables cp, CC, P and Z, we get 


S 2 cp + CjCt + c 2 p = 0 


s - 


K,V 

V 


s + 


V 




= - s a 


(a c + a]_S)cp + b 0 a -(3=0 


w 


( 201 ) 

( 202 ) 

(203) 


or, written in matrix form, 


S - 

b V 


(a 0 + a x S) 


S + 


L 

"v 


K- 


-1 


9 


a 


Sa 


w 


(204) 
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(205) 


A 




The characteristic equation for this system is the determinant 
of the [A] matrix. For this set of equations to have other than trivial 
solutions, the determinant of coefficients must vanish, or 


A = 


K- 


(a Q + a ) 


S + 


K- 

V 


= 0 . (206) 


From equation (206) the characteristic equation is 


S 3 + + c 2 a x ) S 2 + 


C x + 02(^0 + ho) + Aa 1 


+ 


Aao _ Cp k 3 , bp _ 


V 


V 


= 0, 


( 207 ) 


where 


A = 


Cp K p ~ C-| K 5 
V 


(208 


The form of equation (207) is 


B 3 S 3 + B 2 S 2 + B X S + B Q 


= 0 . 


(209) 
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The Routhian array for equation (209) is 


S 3 


B 3 


Bi 


S 2 

S 1 


Bg B 0 

BpBi - B 3 Bq 

b 2 


s° 


Bo 


and, for absolute stability, the coefficients of the first column must 
all be positive. Since B 3 > 0 in our case, the remaining terms are 


B 2 


*2 

V 


K 3 b Q 
V 


i C 2 a 2 


> 0. 


(210) 


BpB, - B 3 Bq = v f C P K 1 + C P K 2 + CxK-A ( 

B 2 (K-l + K 2 ) y + l v J°\ V 


- A a 


+ ( y + C-jCp ) + C 2 3^3 1 9" ( C 2 + 


K A \ 

x ^ 2 j a i -r 0 2 a D a i + ( C 2 + y ) 3 ibo 


( 211 ) 


+ £#2 aobo + C 2 ^ b 2 >0 


B 0 ~ Aa 0 - 



V 


> o. 


( 212 ) 


If at this point we solve for the values of B 2 , (B 2 Bi - B 3 B Q )B 2 1 , 
and B 0 , and find that B 2 < 0 and all other elements of the first column 
are positive, then there is one change in sign between elements B 3 and B^ 
and between elements B 2 and (B 2 B X - B 3 Bo)B£ 1 . Since there are two changes 
in sign, the system has two roots in the right half of the S-plane, and 
the system is unstable. Similar reasoning can be applied for any com- 
bination of signs for the elements of the first column. 
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In the above case, we were able to determine the stability of 
the system for a particular value of gains a 0 , a ± and b Q and system 
constants Cj_, Ki and V. A better picture of the system stability can 
be shown by plotting the stability boundaries in the gain plane (a 0 , b 0 ) , 
(a Q , a-t) or (a x , b D ) . The stability boundaries correspond to the roots 
of S for which the real part is zero. In this case, equations (210), 
(211) and (212) are equated to zero and the equations of the stability 
boundaries are 


b 


K "3 
V ' V 


_l_ — ■ + C^a i = 0 


(213) 


(Ki + K 2 ) -f + 


Sz * , / c^±c^j_c^'y o+ (c^ _ ^ 


\ 


V 


/ 


+ ( ^ + C x C 2 j ai + Cl a Q a ! + (c§ + ]a x b Q 


(214) 


K 3 Cp , in 2 i K^Cp i 2 _ n 
+ a 0 b 0 ' r ^2 As i + y bg - 0 


, CpKn . _ . , Cj 

Aso - b 0 - - b Q + ^ 




0. (215) 


Equations (213) and (214) represent the dynamic stability 
boundary and correspond to the elements of the first column of the 
Routhian array exclusive of the B 0 element. Equation (215) represents 
a stability boundary which can be recognized as the "drift minimum" con- 
trol gain relation. This relation is referred to as the "static" stability 
boundary, or simply the trim condition for stability. 

In order to plot the stability boundaries in the (a Q , b 0 ) plane, 
equations (213) through (215) are written in the following form: 

b D = k iJ (B 2 = 0) (216) 

bo 2 + a obo + j" 3 b Q + a ° + = 0; (B 2 B X - B 3 B D ) B 2 -1 = 0 (217) 
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b Q _ k 6 a Q + k 7 , (B 0 - 0) . 


(218) 


Typical values of the constants for a large booster at the time 
of maximum dynamic pressure are 


C 2 = - .387 
C 2 = 1 . 343 
V = 508 


K x = 20.81 
K 2 = 5.134 

K 3 = 17.66 


and therefore, 


ki 


k 2 


k3 


k 4 


ks 


k6 


k 7 


Cpa i + K JV 
K 3 /V 


K-^Cp 

V 



= -38.66 


= .04671 


= 1.8595 


= 1.8231 


= -.4142 


= .4911 


= +.288 


The stability boundaries, using the above constants and equa- 
tions (216) through (218) for a constant gain of a 4 = 1, are shown in 
Figure 27. We see that the B 2 = 0 stability boundary is a straight line 
with a slope m = 0 and b D intercept 


k i 


f Gpa i + Kp/V \ 
V K 3/V ) • 
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FIGURE 27. 


STABILITY BOUNDARIES IN THE (a , b ) GAIN PLANE USING 
ROUTH'S STABILITY CRITERION ° 
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The B 0 


0 stability boundary is also a straight line with a slope 


K, 


m 


K- 


Ki 


and b c intercept 



The stability boundary for - B-Bo)^^ 1 = 0 is a hyperbola as shown. 

To find the region of stability, we now apply the inequalities given by 
equations (210), (211) and (212). Equation (210) is satisfied for all 
points above the B 2 = 0 boundary, equation (211) is satisfied for all 
points below the B 0 = 0 boundary, and equation (212) is satisfied for 
all points above the upper half of the hyperbola and below the lower 
half of the hyperbola. Therefore, the three inequalities are satisfied 
simultaneously only in the area labeled "stable region." 

Any set of gains (a 0 , bo) for a x = 1 selected from the stable 
region in Figure 27 and used in the given control equation 


(3 = a Q cp + a x q) + b Q a 


would insure that the resulting system was stable. 

From the above discussion, we find that Routh' s stability cri- 
terion provides information about absolute stability of the system and 
very little information about the relative stability. Selecting gains 
from the stable region would not provide information about speed* of 
response, transient decay times, overshoot, etc. 


B. Hurwitz* s Stability Criterion 

The Hurwitz stability criterion is similar to the Routh 
stability criterion in that it provides information about absolute 
stability without actually solving for the roots of the characteristic 
equation. The primary differences between these two criterion are the 
method of expressing the coefficients in the Routhian array and the 
principal minor determinants in Hurwitz 1 s criterion. 
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The roots of the characteristic equation determine the form of 
the transient response and can be obtained from the closed-loop control 
ratio 


M 

R(s) 


k k — l 

S + . . + A-^S 4- A 0 

B n S n + B S n_1 + ... + BnS + B 0 
n- 1 xu 


N(s) 

D(s)’ 


(219) 


where in this case the control ratio is expressed as the ratio of two 
polynomials. The characteristic equation for the system is 

D (s ) = B s n + B S n_1 + ... + B-.S + B 0 = 0. (220) 

n n - 1 -l u 

In a physical system, the coefficients B n , B n _ l3 . .., B 0 are 
usually real numbers. A polynomial equation with real coefficients can 
have either real or complex roots. Any root of D(s) = 0 with a negative 
real part is a stable root and the transient solutions will be of the 
form 


K. e 


-a.t 

x 


or 


-a . t 

K. e ^ (sin u>.t + 0 ) 

J J 

which decay with time* For the system to be stable, all roots of the 
characteristic equation must contain negative real parts. 

The Hurwitz stability criterion establishes the conditions 
under which all the roots of D(s) = 0 have negative real parts without 
actually solving for the roots of the polynomial equation. 
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The n tB order Hurwitz determinant is constructed from the 
coefficients of equation (220) as follows: 


B B 0 0 

n-i n 


0 0 


^ = 


B 


n-3 



B 

n-i 


B 

n 


0 0 

« • 

Bn- B nMl 

0 B o 


( 221 ) 


The principal minor determinants Ai of the Hurwitz determinant 
A n are defined as 


Ai 


Ao 


A 3 


= B 

B 

B 


n- 1 



n-i 

B n 


n-3 

B n -2 


n- 1 


0 

n-3 

B n-2 

B n- l 

n-3 

B n-4 

B n-3 


( 222 ) 


(223) 


(224) 


The Hurwitz criterion then can be stated as follows: The 

necessary and sufficient conditions that all roots of the polynomial 
D(s) = 0 have negative real parts are that B n > 0, A a > 0, A c > > 0, 

An > 0* 

C . Root Locus 


The root locus method was first devised by W. R. Evans in 
1948. Since that time, it has become a very useful method of analysis 
and synthesis of control systems. 

The root locus is a plot of the roots of the characteristic 
equation as a function of gain in the complex S-plane. The poles of 
the transient response mode C(s)/R(s) are related to the zeros and poles 
of the open-loop transfer function B(s)/E(s) and the gain K. The rela- 
tionship is shown below, where G(s) and H(s) are expressed as the ratio 
of two polynomial equations. 
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(225) 


G(s) 


N xi£i 

Di(s) 


H(s) 


Np(s) 

D^Cs) 


(226) 


E(s) 


G(s) H(s) 


Nj_Jk 

D i D2 


(227) 


C(s) 

R(s) 


G (s) = Nj/Dj 

1 + G (s) H(s) Nj N ? • 

®1 ^2 


(228) 


The roots of N x N 2 = 0 in equation (227) are called the zeros 
of the open-loop transfer function G(s) H(s) and the roots of D X D 2 = 0 
are called the poles of the open-loop transfer function G(s) H(s). 
Similarly, the roots of 1 + N 1 N 2 /d 1 D 2 = 0 are called the poles of the 
closed-loop transfer function C(s)/R(s) and determine the operating 
root locations corresponding to the system time constants. The numerator 
N 1 /D 1 of the control ratio C(s)/R(s) merely modifies the coefficients of 
the transient components. 

The characteristic equation of the system is from equation (228) 


1 + G (s ) H(s) = 0 (229) 

or 

G(s ) H(s ) = - 1. (230) 


The location of the roots of the characteristic equation vary 
as the gain K is varied. The plot of these roots as a function of the 
gain K, where K varies from zero to infinity, is called the root locus . 

To establish two important criteria that are useful in plotting 
the root locus, it is convenient to express the transfer function G(s) 
H(s) as the product of linear factors .in the following form: 


G(s) H(s) 


K(S - Z z )(S - Z 2 )(S - Z 3 ) ... (S - Z m ) 

— — n + N ^ m. 

s (s - P X )(S - P 2 )(S - P 3 ) ... (S - P n ) 

( 231 ) 
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In the factored form shown above, the Z T s denote zeros and P's 
denote poles, and pairs of these quantities may be complex conjugates. 
The quantity K is the gain factor of G(s) H(s). For most physical 
systems, the order of S in the denominator is n + N, and is greater than 
or equal to the order of the numerator m. 

Each linear factor of equation (231) can be represented by a 
magnitude and phase angle and is therefore called a phasor. The S-plane 
representation of the phasors, shown in Figure 28, is seen to originate 
on the poles and zeros, and terminates at the point S. 



FIGURE 28, PHASOR REPRESENTATION OF LINEAR FACTORS OF 

G(.) H(s) - Ktt-, 

S (S - PJ (S - P 2 ) (s - P 3 ) 
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The magnitude of a particular phasor, say (S — Pj) where 
S =-a + jw and P x =-ap 1 4- jop.^ will be 


r 


P 


l 



= + N/ r (-a + a ) 2 + 
r 1 



(232) 


and the corresponding phase angle will be 


(J = /S^i. ■ (233) 

p r 

The linear factors in rectangular form can then be related to the polar 
form as follows: 



(234) 
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Upon substituting equation (234) into equation (231), the 
entire G(s) H(s) transfer function can be written in polar form. 


G(s) H(s) 


J n J 7 i^7 

K(r e Zl )(r„ e 2 ) . . . (r e ra ) 


LLl 






. N jN,i P 0 x , Pps 

(r n e °) (r_ e x ) . . . (r_ e n ) 


(235) 




+ 

1 


+'/ 7 )-(N 



+ . . .+/- p ) 

l *n 


= K r e 




where 


r 



(236) 


and 


^ = ( 0 7 + 


+ 


) - (N0 


+ 


+ 


+ 


(237) 


The phase angles ( 0 ) are measured positive counterclockwise from a line 
parallel to the o-axis. 

From equations (230) and (235) then 

G(s) H(s) = K r e^ = - 1, (238) 


where 



cos \|r + j sin >)/. 
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The magnitude of G(s) H(s) is 


|G(s) H(s) j = |K r (cos f + j sin ty) | = | K r j = 1 . 

The gain factor K may be positive or negative. For K > 0, 

G(s) H(s) = K r (cos i|r + j sin \|/) = -1 
= cos \|/ = -1 , 


(239) 


(240) 


since K r = 1 and j = 0 for G(s) H(s) to be a real number. From equa- 
tion (240) the phase angle of G(s) H(s) becomes odd multiples of it. 


t = ± it, ± 3jt, ± 5it, ... 

= (2k + 1)*, k = 0 ± 1, ± 2, ... 


(241) 


For K < 0, 

G(s) H(s) = - Kr (cos \|r + j sin \|/) = -1 


= - cos \]/ = -1, 


(242) 


and the phase angle becomes even multiples of it. 


\J/=0, ± 2 jt, ± 4 tt, ... 

=2krt, k=0, ±1, ±2, ... 


(243) 


In summary, the two criteria that are required in plotting 
the root loci of a system are the phase criterion and magnitude criterion. 

Phase Criterion 

m n 

- ( N0 ?o + I 

i=l i=l 

= 2k it K < 0 

where 


T 


= (2k + 1) it 


K > 0 


(244) 


k = 0, ± 1 


4 * 2 

— ^ 5 • • 
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Magnitude Criterion 


m 

K n r 
i=l 


r 


Po 


n 

n 

i=l 


= 1 , 


(245) 


where II denotes the product of the r's. 

The root locus is plotted by finding all points in the S-plane 
which satisfy the phase criterion. After the locus is completely plotted, 
the magnitude criterion is used to scale it in terms of the values of 
gain K' that correspond to particular roots along the locus. The follow- 
ing procedure is helpful in constructing the root locus: 

1. Obtain the open-loop transfer function and put it 
in the form 


K(S - Zj.) (S - Z 2 ) ... (S - Z ) 

G (s) H(s) = — — ( 246 > 

S (S - PiXS - P 2 ) ... (S - P n ) 

and plot the zerosZ^ and poles in the S-plane 
using the same scale for both the real axis and 
imaginary axis. Continuous curves start at each 
pole of G(s) H(s) for K = 0 and terminate on the 
zeros of G(s) H(s) for K = <». 

2. Draw the loci along those sections of the real axis 
which are to the left of an odd number of real poles 
and zeros when K is positive. When K is negative, the 
locus lies to the left of an even number of real poles 
and zeros. This rule is a result of applying the phase 
angle criterion to a point S as it is moved along the 
cr-axis. The net phase angle contribution of phasors 
originating from complex conjugate poles and zeros to 
the point S on the a-axis is zero and may be neglected. 
For the particular pole zero configuration shown below, 
the real axis portion of the root locus is shown by 
application of this rule. 
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3. The real axis intercepts of the asymptotes are 

located at the centroid of the poles and zeros or 


centeroid 


= 


N+n 

V 

L 

i=l 


m 

v 

L 

i=l 


Z. 

1 


(N+n) - (m) 


(247) 


The number of asymptotes can be found by inspection 
of G(s) H(s) and is equal to the number of poles minus 
the number of zeros or 


number of asymptotes = (N+n) - (m) . 


(248) 


The angle at which the asymptotes intercept the 
a-axis is found from 


0 . 


asymptote 


(2K + 1 ) rc 
(N + n) - (m) 


(249) 


For example, let 


G(s) H(s) 


K(S - 1) 

s^cs + 1) (S + 2 - j2)(S + 2 + j2) • 
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4. The breakaway point between two poles or zeros located 
on the a-axis is shown below in Figure 30. The loci 



FIGURE 30. BREAKAWAY POINTS 


originate at the poles for K = 0 and coalesce at the 
breakaway point a^ and enter the complex region. For 
two zeros located on the a-axis, the locus enters the 
real axis at the breakaway point and terminates on the 
zero at K = oo. in both cases the loci are perpendicular 
to the real axis at the breakaway point. The breakaway 
point can be found by assuming values of a^ between the 
poles or zeros until the following equality is satisfied: 



(250) 


where 

P^g = ith pole to the left of trial point 
P r ^ = i th pole to the right of trial point 
= ith 2 ero to the left of trial point 
Z r ^ = ith zero to t h e r ight of trial point. 
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5. The angle of departure in which the locus leaves a 
complex pole or enters a complex zero can be deter- 
mined by adding up all the phase angle contributions 
to the pole or zero in question. Subtracting this 
sum from 180° gives the required direction. For 
example, the angle of departure for the locus leaving 
pole P 2 in Figure 31 is found from 

0 7 - (0 P + 0 P + 0 P + 0 P ) = ISO 0 


0 p - 180 ° + (0 +0 + 0 p ) - 0 7 . 

*2 *0 *1 *23 



FIGURE 31. ANGLE OF DEPARTURE 


6. The preceding five steps provide enough information 
to construct only a portion of the root locus and to 
determine the manner in which the root locus will 
behave. The remaining branches of the root locus 
located off the real axis can be found by starting 
at a breakaway point and successively selecting trial 
points immediately above the breakaway point until 
the phase angle criterion is satisfied or starting at 
a pole located off the real axis and selecting trial 
points along the tangent defined by the angle of 
departure until the phase angle criterion is satisfied. 
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Once the phase angle criterion for this point is 
satisfied, a new trial point is chosen and adjusted 
until the locus terminates on a finite zero, becomes 
tangent to an asymptote, or enters a short distance 
into the right half S-plane. Only the upper portion 
of the root locus need be determined since it is 
symmetric with respect to the real axis. The lower 
portion of the root locus can be constructed by 
symmetry. Generally, the important part of the root 
locus lies in finite portion of the S-plane around 
the origin. 

The use of a spirule greatly facilitates the 
search of the roots satisfying the phase angle cri- 
terion. It is essentially a device for measuring 
the net phase angle of the phasors drawn from poles 
or zeros to any point S. 

7. Once the locus is completely drawn, it is calibrated 
in terms of the gain K using the magnitude criterion 
of equation (245). 

To show the relationship between the root locus and the cor- 
responding time solution, consider the following second order dif- 
ferential equation relating the attitude cp to the commanded attitude 

<Pc : 


cp + 2 £w n q> + ug <P “ “f! <P C * 


(251) 


The transfer function of this second order system is 

/ \ UJ 

cp(s) = n 

cp c (s) + 2£wn S + * 


(252) 


and the corresponding transient response to a step input on cp c is 


Cp(t) 


S,t 




Spt 


+ 1 


(253) 
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where the roots obtained from the characteristic equation are. 


s i,2 " ± j' jJ n - 1 o ± juj, 0 < £ < 1 

^ l y 2 £ = 1 


S l,2 - 1 > t > 1 

s i,2 “ C = °* (254) 

A plot of the roots S l5 S 2 in the S-plane and the corresponding 
transient solution of cp for a step input of cp c is shown in Figure 32. 

When the roots S l3 S 2 are located to the left of the jco-axis, 
the system is stable, and the transient part of the response decays to zero. 
When the roots are located on the jio-axis, the system contains no damp- 
ing, and sustained oscillations result. As the roots move to the right 
of the jco-axis, the amplitude of the oscillations increases with time, 
and the system is unstable. The above statements are true regardless of 
the order of the system. In the complex S-plane, horizontal lines repre- 
sent constant damped frequency, co^; vertical lines represent a constant 
rate of decay, a; radial lines through the origin represent constant 
damping ratio, £; and circles concentric about the origin represent 
constant natural frequency, con- 


D. Frequency-Response Method 

The frequency-response metnod for the analysis of feedback con- 
trol systems may be described as the study of system behavior with pure 
sinusoidal inputs. It is always assumed that the sinusoidal input has 
been applied for a long time so that the transient response has decayed 
and the steady-state condition has been reached. If the input to a 
linear system is sinusoidal, then the driven response throughout the 
system must be sinusoidal also and of the same frequency, differing only 
in amplitude and phase. The frequency-response approach is the study 
of the amplitude and phase angle of the response of a component or a 
system as a function of the sinusoidal frequency. 

Using the frequency-response approach, the frequency and damping 
of the system cannot be obtained, but enough information can be obtained 
to indicate whether compensating networks are needed. 
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To obtain the steady-state frequency response let S assume 
values equal to jco. The control ratio as a function of frequency is 


C(.jio) _ G(.jgj) 

R(joo) 1 + G(jco) H(jco) ' 


(253) 


For each value of frequency, equation (253) yields a phasor quantity 
whose magnitude is | C ( jco) /R( jco) | and whose phase angle is the angle 
between C(jco) and R(jw). 

The plots in the frequency domain that are used in the analysis 
of feedback control systems are of two categories. In the first category 
are the plots of the magnitude of the output-to-input ratio and the cor- 
responding phase angle versus frequency. These plots may be made in the 
rectangular or logarithmic coordinates. The standard procedure is to 
plot 20 log | G ( j co) H(jco) j (magnitude in decibels) and the phase angle 
versus log co. This plot is often called the Bode plot. The plots in 
this first category are particularly useful in representing the transfer 
functions of individual components. In the second category are polar 
plots of the output-input ratio, called Nyquist plots, and plots of 20 log 
|G(jco) H(jco) | versus phase, called Nichol 1 s plots. The Nyquist plot and 
Nichol's plot have frequency as a running parameter. These plots are 
generally used for the open-loop response. The polar plots when used for 
Nyquist 1 s stability criterion must be drawn for frequencies from -<» to +°° 
(corresponding to the jcj axis of the S-plane). 

In applying the Nyquist stability test, it is required that: 


1. the system be represented by a set of linear differential 
equations with constant coefficients. 

2. the order of the denominator be equal to or greater than 
the order of the numerator of the open-loop transfer 
function G(s) H(s); that is, the 1 im G(s) H(s) 0 

S - > oo 

or a constant. 


E . The Nyquist Stability Criterion 

The Nyquist stability criterion relates the number of zeros 
and poles of the characteristic equation that lie in the right-half 
S-plane to the polar plot of the open-loop transfer function G(s) H(s). 
It can be stated as follows: 
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Given any open-loop transfer function G(s) H(s) , which 

is the ratio of two polynomials in the variable S, let S = jw. 



The amplitude and phase polar plot of G(ju)) H(jcj), as the 
frequency is varied from -oo to +<», encircles the (-1 4- jO) 
point N times. With counterclockwise positive and clock- 
wise negative. 


N = 



Z 


R 


where P R equals the number of poles of G(s) H(s) in the 
right-half S-plane and Zr equals the number of zeros of 
the characteristic equation in the right-half S-plane. 

Since a stable system can have no zeros of the characteristic equation 
in the right-half S-plane, the net number of rotations N of G(s) H(s) 
about the (-1 + jO) point must be counterclockwise and equal to the 
number of poles that lie in the right-half S-plane. 

If G(s) H(s) experiences a net clockwise rotation about the 
(-1 4- jO) point, Zr > Pr and the system is unstable. If there are zero 
net rotations, then Z^ = P^ and the system may or may not be stable, 
depending on whether Pr = 0 or Pr > 0, respectively. If Pr > 0, the 
number of poles in the right-half S-plane can be determined by applying 
Routh 1 s criterion to (denominator of open-loop transfer function 

in equation (227). 

If Nyquist plot representing a system is complicated, it is 
sometimes difficult to determine if the curve encircles the (-1 + jO) 
point, and, if it does, how many encirclements it makes. To determine 
if a system is stable or unstable, trace the curve in the direction of 
increasing frequency. If the (-1 4- jO) point is always to the left of 
the curve, the system is stable. 

If the open-loop transfer function were to pass through the 
(-1 + jO) point, a condition of sustained oscillations would exist, 
representing a marginally stable system. The closer the open-loop 
frequency response characteristic curve comes to passing through the 
(-1 + jO) point the closer the closed-loop system is to being unstable. 
The degree of stability of the system is measured by two factors - gain 
margin and phase margin. 

Gain margin is the measure of the factor by which the gain of 
a system would have to be increased to make it marginally stable. 

Phase margin , 0p^> is a measure of how much additional phase lag is 
required to make a system marginally stable if the gain is unchanged. 

A minimum of 6 db gain margin and 30 degrees phase margin is considered 
to be acceptable margins for system stability in control system design. 
A minimum of 6 db gain margin and 30 degrees phase margin is considered 
to be acceptable margins for system stability in control system design. 
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On the Nyquist plot, the gain margin is the reciprocal of the 
value of G(jco) H(jio) at phase crossover (the point at which the frequency 
response curve crosses the -180° axis). The phase margin is 180° plus 
the phase angle at gain crossover (the point on the plot at which the 
magnitude is unity). 

Using the Nichol's plot (log magnitude versus phase angle), a 
system is stable if the 0-db, -180° point is always to the right when 
tracing the curve in the direction of increasing frequency. Gain 
margin and phase margin can be determined quite readily from the 
Nichol's plot. The gain margin (in decibels) is measured where the 
curve crosses the -180° axis. The phase margin is 180° plus the phase 
angle at the point where the curve crosses the 0-db axis. Quite 
often the magnitude of the response (in decibels) is plotted versus 
phase angle plus 180°. For this plot the origin is the 0-db, 0 phase 
margin point instead of the 0-db, -180° point. This plot enables one 
to read phase margins at a glance. 

On the Bode plot the gain margin is measured in decibels on 
the log magnitude curve where the phase-shift curve crosses the -180° 
line. Phase margin is measured on the phase-shift curve at the fre- 
quency where the log -magnitude curve crosses the 0-db line. 

The easiest step in adjusting for a satisfactory system response 
is to adjust the gain. If satisfactory response cannot be achieved by 
gain adjustment alone, then compensating techniques must be used. 

A root locus plot, Bode plots, Nyquist plots, and Nichol’s 
plots have been drawn for a system of equations to illustrate many of 
the things that have been discussed. The following set of equations, 
considering only one bending mode, was used: 


Cp = - G-^OC - Cop 

(254) 

a = 9 

(255) 

cp f = cp + ^Y’Cx^ 

(256) 



$ f = 9 + )/• Wi 

(257) 
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(258) 


T). + 2t. W . Tj . + CO? T) . 

b l 1 *1 1 ‘l 


R'Y.&gp + p Y.(x h ) - I E Y^) 


m . 

l 


6 +2 1 co 3 + co 2 3 = co 2 e 

c ’c c c c c c 


g = a D cp + a x cp - a Q cp c 


P A " T ^ c " ?A> 


3 + 2^ e ^ 3 + <o| 3 = «*! e, 


(259) 

(260) 
(261) 
(262) 


By substituting equations (255), (256), and (257) in the other equations, 
the system equations are reduced to six equations in six unknowns. 


cp “ c x cp — c 2 3 


(263) 


n . + 2 1 . co. n . + (o? rj . 

'i h i 'i i 'i 


R,Y dv p + Ve w • 


i5 


m . 


(264) 


6 +2 1 co 3 + co 2 3 = co 2 e 

K c b c c c c c 


(265) 


e 


a CD 
cr 


+ a 

o 


X Wi + a i(p + a i 



(xcp)^ - a Q cp c 


(266) 


Pa = T( Pc - P A > 


(267) 


3 + 2^ e cog, 3 + 3 = 3 a - 


(268) 


117 



The Laplace transform equations are 


(S 2 + C x )cp = -C 2 £ 


2^ = 


(S^ + 2t. go.S 4- w f)q 
1 i 


R ' Y i (x h> + s£ <Ve W - r E W 


m . 

l 


(S 2 + 2£ go S + cj 2 )P = go 2 e 
c c c c 


€ = (a x S + a Q )cp + 


a i s Y Y^(xcp) + a Q y Y^(xcp) 


T] . - a cp 
o Y c 


(S + T)3 a = Tf3 c 


(s * + 2 ^e “e s + = °| \ 


(269) 

(270) 

(271) 

(272) 

(273) 

(274) 


where 

= control filter damping ratio 
go c = control filter frequency 
T = actuator lag time constant. 

By combining equations (273 and 274), the equations can be reduced to 
five equations in five unknowns. 


(S 2 + C x )9 = - 

(S 2 + 2f . CJ.S + 
^1 1 

(S 2 + 2t CJ S + 
s c c 


c 2 5 



R' Y 




i (x h ) 


+ S 


2 (m E i„ Y 


i (x h> 


m . 

i 




(275) 

(276) 

(277) 


118 



€ 


( a l S + a 0 )$ + 


a iS 2 (x<p) + a Q Yl (xcp) ^ 


a cp 
o Y c 


(S + T) (S 2 + 2^ E 0i E ,S + 



■ H K- 


These equations may be represented by the following block diagram: 


Rigid Body 



FIGURE 33. BLOCK DIAGRAM OF FEEDBACK CONTROL SYSTEM 







where 


G E (S + T)(S a + 2£ EUfe S + ug) 


G c + IF to S + uf 1 
c c 


. .Sz — 

+ cl 

R> W + S 2 ( m E ^ E Y.(x h ) 

L m. (S^ + 2t .cj.S + (Jr) 

1 b i 1 1 


X E W 


A - a L S W- ( V *», [ Y L<y} 


The block diagram shown in Figure 33 can be reduced to the basic single 
loop feedback control system as follows: 

(a) Move the flexible body feedback loop into the rigid 
body feedback loop. 







(b) Combine the elements in the feedback loop. 



The above block diagram in (b) can now be reduced to the following form: 
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Open-Loop Transfer Function 


9 f 

~ = G(s) H(s) 

e 


= G G 
c E 


(a + a,S) G + A G 

o ^ Cp T] 


T to 2 (B + C) 

(s + THs^+CiHs 2 + 2^0^ +to 2 E )(s 2 + 2^1Ts + topes* + 2£ u s + to 2 ) 

(288) 


where 


B = 02(3!$ + a D )(S 2 + 2£ w S + to?) 


C = 


a 0 > Y^) + ai s ') V (x,p 


! ' Y i (x h> 


+ S £ 


Ve W- I E r(x h ) 


m. 


Closed-Loop Transfer Function 


9 

a o9c 


G(s) 

1 + G(s) H(s) 


(289) 


9 

% 


a G 

o c 

1 + G G_, 
c E 


G 

-9. 


[ (a Q + aiS) 


G cp + A G^] 


(290) 


a C 2 Tto 2 to 2 (S 2 + 2f . to. S + to?) 
o K c b l i i 

D i + Do + D 3 
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where 


D x = (S + T)(S 2 + C 1 )(S 2 4- 2 C e w e s + ^)(S 2 + 2 ^c L °c S + ^c^ 52 + 2 ^i W i S + 
D 2 = Taj 2 co 2 C 2 (a iS + a Q ) (S 2 + 2£^iihS + w 2 ) 


D 3 = Toj| (J 2 (S 2 + C x ) 


a nS 


y ; (X.) + a 

1 cp ° 


n<v 


R ' W + s *<Ve W - *e W 


m . 


Characteristic Equation 

1 + G (s ) H(s) = 0 


or 


D i + D 2 + D 3 = 0. 


(291) 


(292) 


The closed-loop transfer function and the characteristic equation 
can also be found by expressing the system equations (equations 275 
through 279) in the following matrix notation: 


A X = a Q U 

where 




I 

0 

1 

p 


0 

n 

II 

1 ^ 

0 

h 


*c 

e 


0 


(293) 


(294) 
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0 


0 


0 


S^fCi 


c 2 


R'Y (x )+S 2 (m £ ■ Y (0-1 Y' (jc ) 

0 \ E E i _h E 1 11 . / -(S^t.w.S+uf) 0 

1 1 


m . 


A = 


0 (S £ +2t oj S+io 2 ) 

b c c c 


- OJ 




( a i s+a o) 


a x S ) Y\ (x. )+a Y'. (x ) 

/_ j 1 9 ° cp J 


-1 


(S+T) (S 2 +2^ eWe S+u|) 


-Tuf; 0 

E 


(295) 

The closed-loop transfer function cp / Cp c may be obtained by setting 


A il 

(296) 

All 


where A' x is the matrix with a Q U replacing the first column (the column 
describing the contribution of cp to the homogeneous equation). 


Co 


R 'Y i (x h ) + S 2 ( m E iE Y i ( x h )- I E Y'i(x h ) 


m. 

L 


- ( S + 2 f . Cl) . S+ GO . ) 0 

^111 


9 


0 (S 2 +2f co S+uT) -co^ 

^c c c c 


l O ( ?C 


a jS ) Y r . (x. )+a 0 Y f . (x ) 

x ' 1 1 cp 


-1 


0 (S+T) (S s +2^ E aj E S+ca|) 


-H 
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(297) 



By determinant expansion of the numerator and denominator, the closed- 
loop transfer function becomes 


Cp Tcoj; (S^ + 2f^utS -I- oj?) 

Cp c D-j_ + D 2 + D 3 


(298) 


which is the same as the closed-loop transfer function obtained from 
the block diagram. 

The characteristic equation is given by + Do + D-, which has the 
same roots as one plus the open-loop transfer function. The open-loop 
transfer function is given by (D 2 + D 3 )/D 1 , which is the same as the 
open-loop transfer function obtained from the block diagram. 

To obtain the root-locus plot, the characteristic equation was 
programmed on a digital computer. The computer calculates the coef- 
ficients of the characteristic equation as the gain a 0 is varied from 
zero to very large values, and then solves for the roots of the result- 
ing polynomial equation. 

The root locus plot for a typical large booster at the time of 
maximum dynamic pressure is shown in Figures 34 and 35. The system is 
initially unstable when the gain a Q = 0 since one root is located at 


+ \ - Cj 


+ 





where in this case the aerodynamic center of pressure location (x C p) is 
forward of the vehicle center of gravity location (x C g) , making this 
root lie in the right-half S-plane. 

This system remains unstable until a Q is increased to .25. The 
system again becomes unstable for a gain a Q > 1.5 due to the pair of 
complex conjugate control filter roots entering the right-half S-plane. 

The poles of the open-loop transfer function could also be obtained 
from the charac ter is tic equation since 


1 + G (s ) H(s) = D x + D 3 + D 3 = 0 


= 1 + 


± P; 

D n 


0. 


( 299 ) 
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Then the open-loop transfer function is 


G(s) H(s) = . (300) 

U 1 

It was shown in Section IX (C) that the roots of D-l = 0 are the 
poles of the open-loop transfer function; therefore, 

Dx = (S+T) (S 2 + CjMS 2 + + wj;) (S 2 + + w 2 ) (S 2 + 2( u S + a) 2 )=0, 

(301) 

and 

S = - T 

S = ± N P C x C x < 0 


s ■ - Ee“e * H" !l - tl t E < 1 


s - - C c « c ± j« c -' 1 - Ce ?c<! 


s = - e.w. + jw. vi - t? . £ < i 

These nine roots are the starting points for the root locus when 
the gain a Q = 0, and are shown in Figures 34 and 35. 

The roots of D 2 + D 3 = 0 in equation (300) are the open loop zeros 
as shown in Section IX (C) . In this case D 2 + D 3 = 0 can be written 
as the following fifth order polynomial: 

B 5 S 5 + B 4 S 4 + B 3 S 3 + B 2 S 2 + BjS + B q = 0, (302) 


where 

B s = -^- L Y' (x.) K 4 

5 a D cp 
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'- Y '<y k * 


= T* 1C 2 - Y' (x ) [C^ + k 5 ] 
a o [ cp 


Bp = 


^ + c 2 - Y' (x ) [C x K4 + K 5 ] 


Bn = 


= 2 t .go. 
1 


Cp -f 


ai 


C 2 wf - Y' (x) K 5 C x 


B 0 = C 2 cj? - Y' (x ) K5C1. 

1 cp ^ ■ L 


K, = 


Ve W ■ y e W 


R 


K. - 


W 


m , 

1 


m. 

l 


For typical values of the constants in the coefficients, 
equation (302) becomes 


S 5 + 3.85 S 4 - 200.45 S 3 + 1462.45 S 2 - 20710 S - 18785 = 0, 


(303) 


and the five roots are 
S * 13 
S = -20 
S = -.85 
S = 2 ± 9i. 

These are the points on which the root locus terminates as the 
gain approaches infinity. They are shown as small circles in Figures 34 
and 35. Since there are nine open-loop poles and five open-loop zeros, 
four of the loci must terminate on zeros located at S = 

To illustrate Routh 1 s criterion, the characteristic equations 
for an a 0 of 1.4 and 1.5 were set up in the Routhian array (Table 4). 

For an a Q of 1.4, all signs in the first column are positive, indicating 
a stable system; for an a G of 1,5, there are two sign changes in the 
first column, indicating two roots in the right-half S-plane and hence, 
an unstable system. 
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FIGURE 35. BLOWN UP VIEW OF ROOT LOCUS SHOWN IN FIGURE 34 IN THE REGION OF 
THE ORIGIN 
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TABLE 4 

ROUTHIAN ARRAY 


a o 

= 1 , 

.4 


+1. 

,000 

+00 

+2 

+ 2. 

,864 

+01 

+4 

+2. 

,632 

+04 

+5 

+3, 

.714 

+04 

+1 

+4, 

.534 

+06 

+1 

+5, 

,887 

+05 

+1 

+7, 

.155 

+07 

+7 

+4, 

.227 

+06 

+2 

+2, 

.599 

+08 


+2, 

.626 

+07 



a o = 

.5 


+1.000 

+00 

+2 

+2.864 

+01 

+4 

+2.632 

+04 

+5 

+3.714 

+04 

+1 

+4.531 

+06 

+1 

+5.691 

+05 

+8 

+8.452 

+07 

+8 

+2.900 

+06 

+3 

-1.075 

+08 


+3.126 

+07 



436 

+03 

+2.722 

347 

+04 

+1.982 

815 

+06 

+1.611 

807 

+06 

+1.745 

00 

+08 

+9.068 

002 

+07 

+2.626 

045 

+08 


626 

+07 



436 +03 +2.720 

347 +04 +1.981 

810 +06 +1.627 

804 +06 +1.694 

506 +08 +1.052 

312 +06 +3.126 

036 +08 
126 +07 


+05 

+6.269 

+06 

+06 

+1.845 

+07 

+08 

+9.254 

+08 

+07 

+2.626 

+07 

+08 



+07 




+05 

+6.312 

+06 

+06 

+1.811 

+07 

+08 

+1.074 

+09 

+07 

+3.126 

+07 

+09 



+07 




+3.322 +07 
+2.626 +07 


+3.861 +07 
+3.126 +07 
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For the Nyquist plots (Figures 36 A-C), Nichol’s plots (Fig- 
ures 37 A-C), and Bode plots (Figures 38 A-C), three different gains 
were used: .2, .7, and 1.5. Gains of .2 and 1.5 represent unstable 

systems; gains of .7 represent a stable system. Gain and phase margins 
are shown on the Bode plots, Nyquist plots, and Nichol’s plots where 
it is possible to indicate these. The system considered is a con- 
ditionally stable system; that is, it can be made unstable by increasing 
or decreasing the gain. Since this is true, there are two gain margins 
to consider - a factor by which the gain is increased and a factor by 
which the gain is decreased. On the Bode plots and Nichol’s plots, a 
positive or negative sign has been used to indicate whether the gain 
should be increased or decreased. 



FIGURE 36A. NYQUIST PLOT, p = . 2cp + .22cp 
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X. ADAPTIVE CONTROL 


As a booster vehicle moves on its trajectory from the dense atmosphere 
at low altitudes to the rare atmosphere at high altitudes, extreme changes 
occur in the aerodynamic parameters and mass distribution. A fixed set 
of gains and compensation time constants may be incapable of providing satis- 
factory performance throughout the trajectory. Instead, a gain scheduling 
procedure is required in the control elements as the surrounding condi- 
tions change. This procedure is not particularly attractive for several 
reasons . 

First, accurate information is required about the aerodynamic and 
structural characteristics of the vehicle and its environment. Second, 
there is the need to measure air data. Third, the establishment of the 
gain scheduling is a long difficult process. 

The attractive features of a self-adaptive control system are that 
it eliminates the need for accurate information about the controlled 
system and it adjusts the control parameters on the basis of the amount 
of deviation existing between the actual and desired output response. 

Various types of self-adaptive control systems have been developed 
in recent years. No attempt is made to discuss the various types of 
self-adaptive control systems, but only a description of one of the most 
prominent systems which is applicable to the control of aerospace vehicles 
is presented. 

The following adaptive control system was developed by Minneapolis- 
Honeywell [3] for a large, highly elastic booster vehicle. The two out- 
standing virtues of this system are its inherent ability to self-adjust 
its control gain to compensate for changes in control moment and its 
tolerance to variations in aerodynamic characteristics and structural 
mode shapes and frequencies . 
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Model Concept 


Adaptive Gain Changer 






A. Model Concept 


The model concept utilizes a simple lag filter (the model) to 
shape the attitude error signal which serves as a command to a relatively 
fast attitude rate inner loop. By making the band width of the attitude 
rate loop at least three times larger than the bandwidth of the model, 
the over-all response between attitude error (attitude rate command) and 
vehicle attitude rate is essentially established by the model. 

The block diagram shown in Figure 40 can be reduced to the 
simplified block diagram as shown in Figure 41. 

The closed-loop transfer function is obtained by solving the 
following system of equations obtained from the block diagram in Fig- 
ure 41 for the ratio cp/cp c . 


= M(s) cp c 

(304) 

e = 9m “ 9f 

(305) 

<p = K c A(s) e 

(306) 

9f = B(s) cp. 

(307) 


The closed-loop transfer function then is 


K A(s) M(s) 

Cp c 

<P ~ 1 + K c A(s) B(s) 


A (s ) M(s) 

-t + A (s ) B (s) 
^c 


(308) 


If the control gain K c in equation (308) is increased to large 
values the transfer function between an attitude change and an attitude 
command approaches a fixed second order characteristic: 


1 im 

K c -> » 


9c 



K 


JL 


T S 
m 


t)C 


T S + 1 
m 


T m s^r 


S + 



s 2 + 


K cp /T m 

— S + % 


m 




(309) 
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FIGURE 40. NORMAL SYSTEM CONFIGURATION 


M(s) A(s) 



FIGURE 41. SIMPLIFIED BLOCK DIAGRAM OF ATTITUDE CONTROL SYSTEM 


139 










Here it can be 


seen that the desired frequency 



and damping 


i 


m 


1 

2 



cp m 


are controlled exclusively by the model time lag and attitude gain 
K<p if the control gain K c is maintained at a sufficiently high value. 

The action of the control system in forcing the vehicle toward 
the desired frequency and damping is illustrated by considering the locus 
of the closed-loop poles of the system as a function of changes in gain 
K c as shown in Figure 42, It is evident from this plot that increasing 
control gain forces the unstable vehicle poles toward the pair of zeros 
established by the combination of rate feedback plus the fixed model 
time constant and fixed attitude gain. These zeros are established at 
the desired natural frequency and damping ratio by choosing appropriate 
values of Tm and K^p for K c oo. For this case the frequency and damping 
was chosen to be 


Frequency: 1.87 



Damp ing : 


.68 


1 pH 

2 -h m K 0 


then, 


K0 = 1.4 


TL = 
m 


.4. 
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The open-loop transfer function is 


K c A(s) B (s) 


K C 2 [S + (C 2 K 2 - C^)] (S 2 + ^ S + ^) 

c LpV 1 m x m 

(S 3 + S 2 + C X S - + T_) 


(310) 


Using the values of *C l3 C 2 , K 2 , K 3 and V given in Section IX-A and 

the values of K 0 and T m above, equation (310) becomes 

1.343 K c (S + .0219) (S + 1.25 + 1.39j)(S + 1.25 - 1 . 39 j ) 
K c A(s) B(s) = (S + .647) (S - .596) (S - .042)(S + 2.5) 

(311) 

The zeros are contained in the numerator and the poles are con- 
tained in the denominator of. equation (311). These are plotted in 
Figure 42. The arrows indicate the direction of movement of the closed- 
loop poles as the gain K_, is increased. 


jco 



FIGURE 42. ROOT LOCUS OF SIMPLIFIED ATTITUDE CONTROL SYSTEM 
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If the control gain K c can be made sufficiently large, the basic 
system is highly tolerant to changes in vehicle characteristics as can 
be seen from equation (309) where the rigid body transfer function com- 
pletely drops out. It becomes increasingly difficult though to achieve 
the necessary high inner loop gain as additional control lags (low fre- 
quency structural bending modes) are introduced in the control loop, 
unless some means of controlling the loop gain is employed. 

B. Adaptive Gain Control 

If the vehicle aerodynamic and structural data are known 
accurately and if there is adequate separation between the desired 
rigid body frequency and the lowest structural mode, adequate control 
can be achieved by scheduling control gain in an open loop manner, e.g., 
changing gain as a function of time. Basic to this approach is the 
requirement for adequate rigid-body and structural mode gain margins. 
Current trends in the design of booster configurations make this 
increasingly difficult to accomplish. 

With adaptive gain control, specified gain margins are 
unnecessary because the gain level established in the system is based 
on measured inflight characteristics rather than predicted values. 

The adaptive gain changer operates to maintain the inner loop gain 
either at its maximum stable value or at some fixed percentage depend- 
ing on the mechanization employed. 

The rate loop stability is discussed first to establish the 
basic premise on which the gain changer concepts are founded. For 
simplicity, only the rigid body will be considered and the results 
are applicable in principle even when all system dynamics are included. 

The simplified rate loop block diagram with unity feedback is 
shown below. 



FIGURE 43. SIMPLIFIED BLOCK DIAGRAM OF RATE LOOP 
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(312) 


2 = 
e 




S + K 


-CpS 


+ dS + c 


ff = 


e 


- C 2 (K K S + K. K ) 
^ c p 1 c 

+ dS + C ± 


G (s ) . 


The closed-loop transfer function is obtained as follows: 


(313) 



? = 9 - 9 m 

(314) 

cp = G (s) e = G (s) (9 - «p m ) 

(315) 

cp (1 - G(s) ) = - G (s ) <p m 

(316) 

£_ - G ( s > 

7 1 - G(s) 

(317) 

%i 


Substituting the open loop 

transfer function from equation 


(313) into equation (317) and rearranging we get the closed-loop 
transfer function of the rate loop. 


SL = 


(K K S + K.K ) 
C 1 c 


+ (d + CpK K )S + ( C *i + CoK.K ) 

*— r* n — ^ t C 


c p 


(318) 
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The characteristic equation for equation (3 1 8) is 


S 2 + (d + C^K K ) S + (Cl + C-K.K ) = 0. (319) 

C P y X ^ L C 

For large values of K c C 2 and selected Kp and K x , the terms 

N 1 

d = — — = aerodynamic damping 


C x = aerodynamic moment coefficient 

may be considered negligible. The characteristic equation then becomes 

S 2 + Co K K S + K K. = 0 . (320) 

^ c p ^ c 1 

Maintenance of K c at a value where the approximation of equation (320) 
is valid is essentially the objective established in the development of 
the gain changer concept. 

The form of the characteristic is 


S 2 + 2£u>S + OJ 2 = o, 


(321) 


and by comparison of terms, it is^ seen that the bandwidth of the rate 
loop is given by 


gj = Co K.K 
r ic 


and the damping by 


5 


C-3 K K K 

1 c P = _R 

2u 2 



(322) 


(323) 
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where 


C 2 = control moment coefficient 
K c = forward loop gain 
Kp = forward proportional gain 
= forward integral gain. 

Except for the term C 2 , the rate loop dynamics dominate all inherent 
characteristics of the vehicle. By referring to equation (320), it is 
evident that completely uniform characteristics are maintained by an 
inverse variation of loop gain K c with C 2 . 

Loop Gain K c Operating at Critical Value 


Operating the loop gain K c at its critical value forces the 
system into oscillatory motion (i.e., the operating poles move to the 
jojaxis). The oscillation develops into a limit cycle having a con- 
stant small amplitude at all flight conditions. From equation (322) 
it is seen that maintaining K c at a maximum value provides the maximum 
rate loop bandwidth and maximum rigid body stability margin. The limit 
cycle frequency is determined by the frequency of the least stable root 
locus crossing into the right half plane. On most large vehicles the 
first bending mode will establish the practical limit to the bandwidth 
of the attitude rate loop. The limit cycle to operate the gain changer 
will then be close to the first bending mode frequency. 

For a given sensor location and control dynamics, a gain cor- 
ridor will be available similar to that shown in Figure 44. The width 
of this corridor will depend on the compensating networks and the loca- 
tion of sensors. The upper gain limit, which represents the high fre- 
quency stability boundary, will be primarily determined by the attitude 
rate loop dynamics. The lower gain limit, which represents the low 
frequency stability boundary, will be determined by both the inner loop 
(attitude rate) and outer loop (attitude and lateral acceleration) dynamics. 

As the bending mode frequencies decrease, the upper gain limit 
moves down; as the degree of static instability increases, the lower gain 
limit moves up. The use of the automatic gain changer will keep the sys- 
tem operating along the upper gain boundary and will provide maximum 
rigid body stability margin. The operation of the gain changer is 
described below with reference to the block diagram in Figure 45 and 
root locus in Figure 46, 

Input to the unit is applied to two filters from the signal 
representing engine deflection. The down logic filter is a band pass 
which passes frequencies near the chosen limit cycle frequency. The up 
logic filter passes frequencies near the rigid body frequency. 
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FIGURE 44. TYPICAL AVAILABLE CONTROL GAIN ALONG TRAJECTORY 
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FIGURE 45. GAIN CHANGER BLOCK DIAGRAM 
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FIGURE 46. SIMPLIFIED ROOT LOCUS 









The output of the down logic filter is rectified and subtracted 
from a constant voltage called the set point voltage E^. Any resultant 
voltage Eg > 0 is integrated and applied to a multiplier which is in 
series with the missile control loop so that the loop gain will rise. 

The loop gain will continue to rise until the output of the down logic 
filter cancels the set point voltage. Thus the set point voltage 
establishes the stable limit cycle. 

If sustained disturbances are present which contain large 
amplitude frequencies near the down logic band-pass frequency, then 
it is possible that the gain will be driven below allowable limits 
dictated by rigid body stability margins. In this case any resulting 
low-frequency rigid body oscillation is sensed by the up-logic filter 
and its output voltage is summed with the set point voltage to drive the 
gain up. 


A bracketing of gain is thus achieved so that a compromise is 
obtained between the allowed amplitudes of high and low frequencies for 
all disturbance conditions. 

The rate at which the gain is changed can be controlled by 
adjustment of the proportional plus integral device in conjunction with 
the limiter. These limits serve to minimize the effects of large 
actuator excursions on the gain changer. The filter sets the static 
gain of the gain changer and reduces transmission of large transient 
voltages and rectifier ripple. 

C. Gyro Blender 

The gyro blender was developed so that a favorable total rate- 
gyro signal for the first bending mode may be obtained at all times of 
flight with minimum positioning of the gyros. One gyro is located 
forward and one aft of the first bending mode antinode. The outputs 
of two rate gyros are automatically blended to establish a desired 
magnitude and phase relationship of the first mode pickup. The blender 
is designed so that it will not affect the rigid body output but will 
give a positive, negative, or zero signal for the first bending mode 
depending on the system requirements. 

Figure 47 is a block diagram of the gyro blender. The blend- 
ing of the gyro signals is accomplished by a proportional servo used 
to position two potentiometers. 
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FIGURE 47. GYRO BLENDER BLOCK DIAGRAM 

The forward and aft gyro signals are 
3 

= *R + Z 

i=l 
3 

‘ + X 

i=l 

where 

<Pp, = forward and aft rate signals 
Cp^ = angular rate of change of missile rigid body 
tnode shape 

Tj^ = i t ^ 1 mode generalized amplitude function. 


(*F (324) 

(Y; (325) 
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Depending on the position of the blender input potentiometers, the 
summed output sent to the autopilot is 


= (1 - K) 9 + K cp. 


'A T F 
3 


= (1 - K) 


iP R + Z, (Y a 

i=l 


+ K 


(326) 


+ > <n ti) 


F i 


i=l 




(1 - K) Y' + K Y' 
lA iF 


(327) 


i=l 

3 


= + Z W '°i 


(328) 


i=l 


where 


A' =n: [t (i. K) 


(329) 


and K is the blender potentiometer position expressed as some fraction 
of unity. 


Positioning of the blender potentiometer is accomplished by 
rectification of two band-pass filters which peak at or near the first 
bending mode frequency. Since the inputs to the filters are the gyro 
outputs modified only by attenuator K x the blender potentiometers will 
be driven to a position where the input to the integrator is zero. 
Since the filters pass only the first bending mode frequency and both 
signals at the summing junction of the filters are equal, we can write 
the expression below: 

|K X K Y’ 1F | = |1 - K| |Y' 1a | , (330) 
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where K x is the fraction of voltage passed by the attenuator, assuming 
the blender filters pass only the first mode frequency. An expression 
can be derived for A^ which is the first bending mode slope sent to 
the flight control system. By eliminating K from equation (329) using 
equation (330), the bending mode blended slope is as follows: 


Ai = 


(1 - K ) Y* Y* A 
v r if iA 


K Y' 

1 lF 


+ |Y' 


lA 1 


(331) 


Since in this case the attenuator K x has been placed in the 
forward gyro input to the blender, the blender potentiometers will be 
positioned so that a greater portion of the forward gyro input relative 
to the aft gyro output will be sent to the control system. Thus, by 
placing K x in either one or the other of the aft or forward gyro inputs 
to the blender, the magnitude and sign of Aj_ can be controlled. 

The magnitude of A^ affects the separation of the first mode 
poles and zeros, while the sign of Aj determines the position of the 
zero above or below the poles in frequency. If the potentiometer K x 
is placed in the forward gyro input to the blender, the first bending 
mode zero frequencies will be greater than the pole frequencies. The 
exact position of the first bending poles and zeros in an actual mis- 
sile system will be influenced by coupling effects. However, with 
the stipulation that the blender filter circuits have sufficiently 
attenuated all other frequencies, the blender can be made to position 
the first mode poles and zeros to give the desired orientation. 

The limiter in conjunction with the gain K 2 will limit the 
rate at which the blender can be changed. Thus, short duration fre- 
quencies other than the first bending frequency will have little 
effect on the blender. 
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APPENDIX A 


ACCELEROMETER EQUATION 
(Reference 9) 



The basic elements of a linear accelerometer are shown above in 
Figure 48. The accelerometer is attached to the vehicle structure and 
senses the resultant acceleration Zf . The axis of the accelerometer is 
oriented normal to the vehicle longitudinal axis and measures the 
accelerations due to lateral motion. 


- Z - (x 


eg 


" *a> 


9 + 


\ w» 


(Al) 


and due to the lateral 
or. 


component of the longitudinal acceleration 




l 


(A2) 
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The resulting acceleration of the accelerometer frame is 


z _ = 


+ z P = Z 



x a )$ + 



9 


Y . (x ) . 

1 a ‘i 



V (x a ) 



(A3) 


The accelerometer equation of motion becomes 


m 2- 4* C z + K z = - m z 
a a aa aa af 


or 


C K 

.. . a • . a 

z + — z + — z 
a m a m a 
a a 


2 

f 


where 


C 

t = — 
b a 2m oo 
a a 



(A4) 


(A5) 


The accelerometer's frequency is designed to be very high by 
virtue of a small mass and a very stiff spring. Accordingly, there 
would be very little displacement if it were used to detect low fre- 
quency -oscillations in the missile airframe. It is convenient to intro- 
duce the scale factor 
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to measure frequencies which are small in comparison with the acceler- 
ometer natural frequency. Introducing equation (A6) into equation (A5), 
the accelerometer reading becomes 


A 2t A 

-§ + —2 — 2. + A = z 

oo to a f 

a a 


or 




A 

a 


- z 


f 


and the accelerometer transfer function is 


(A7 ) 


(A8) 


A 

a 




(A9) 


From equation (A9) for to « ooa, it is seen that the acceleration 
sensed by the accelerometer is approximately equal to the input accelera 
tion or 


f 


A . 

a 


(A10) 


The block diagram relating the input acceleration to the sensed accelera 
tion is shown below. 



FIGURE 49. ACCELEROMETER BLOCK DIAGRAM 
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APPENDIX B 


RATE GYRO 
(Reference 1) 



FIGURE 50. RATE GYRO 


The figure above shows the basic elements of a rate gyro. In 
many control systems, it is sometimes necessary to measure the angular 
rates of the vehicle. For this purpose the rate gyro is used. 

The inner gimbal supporting the gyro is restrained by a spring 
which permits an angular displacement 0 O . The outer gimbal is rigidly 
attached to the vehicle structure. As the vehicle rotates about the 
input axis at an angular rate cp f , the rate of change of the angular 
momentum vector H is H<^£ which provides a moment M about the output 
axis. The equation of motion of the rate gyro then is 


J0 O + C0 O + K0 O = H q> f 


or 


(S2 + J S + 00 = J $£ 


(Bl) 


(B2) 
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where 


i 


g 


c 


2J u) 

g 


CJ 

g 



(B3) 


(B4) 


The transfer function between the angular displacement and the 
angular rate is then 



H/J 

S 2 + 2f co S + co 
g g g 


(B5) 


In most cases the angular displacement Q (J is related to an output 
voltage E Q , using a potentiometer, where 


Eo ” K v Bo* 

The block diagram for the rate gyro is shown below. 


(B6) 



FIGURE 51. RATE GYRO BLOCK DIAGRAM 


E q = output voltage 
C = viscous damping 
J = gimbal inertia 
K = spring constant 


Ky = pickoff sensitivity 

0 O - output angular displacement 


input angular rate q) + Z f]^ Yl (x ) 


H = angular momentum. 
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APPENDIX C 


RATE INTEGRATING GYRO [1 ] 



The rate integrating gyro is used to establish an attitude reference 
for the control system. The vehicle’s course or attitude can be changed 
by programming flight control commands into the command torques from a 
computer. The deviation of the vehicle longitudinal axis from the atti- 
tude reference provides the necessary error signal to command a deflec- 
tion of the thrust chamber. 

% 

The integrating gyro as shown in Figure 52 is designed so that any 
turn of the gimbal works against the torque of the viscous restraint, 

C0 O . When the integrating gyro senses an angular rate ^ about its input 
axis, the gimbal has to overcome the torque of the viscous restraint in 
order to precess. 

Unlike the spring restrained rate gyro, the gimbal will continue to 
prfecess as long as there is an input. The resultant effect of a distur- 
bance is to produce an angular rotation about the output axis that is 
proportional to the integral of the input axis rate. The input distur- 
bances come from two sources: command torcues and input axis rates. 

The equation of motion is 

(JS 2 + CS) e o = Hq> f + K t E t , (Cl) 
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and the respective transfer functions are 


Command Torques : 



i K K 

i p 


S JS + c 


(C2) 


Input Axis Rates : 


E 

o 



-I HK 

i E 

S JS 4* C 


where 


(C3) 


E 0 “ K p 9o- 


(C4) 


The total transfer function is 

(H9 £ + K t Ei). (C5) 

S + I 


The block diagram of the rate integrating gyro is shown in 
Figure 53. 
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Torquer 



E q = output voltage 
C = viscous friction 


H = angular momentum 

J = gimbal inertia 

Kp = pickoff sensitivity 

K t = command torques gain 

0 O = output angular displacement 


<? f 


input axis rate. 



FIGURE 53. BLOCK DIAGRAM OF A RATE INTEGRATING GYRO 
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APPENDIX D 

AERODYNAMIC VANE TYPE SENSOR 
(Reference 9) 



FIGURE 54. ANGLE OF ATTACK METER 


The aerodynamic vane sensor is used to measure the angle of attack 
between the vehicle longitudinal axis and the direction of the relative 
velocity vector. The moveable vane type sensor is located at the nose 
of the vehicle as shown in Figure 54. 

From Figure 54 the following relations are obtained: 


q) a = cp + ^ Y [( x a> - °v ' (D1) 

i 
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a = ql , - a 

a i v 


(D2) 


a = a 

T 


+ I W \ - a r 


(D3) 


The angle aj_ is an induced angle, due to the bending and attitude 
motion, given by 


“L 


X Y i (x cP \ + (x cg ' x cc ) $ 


(D4) 


and the angles are defined as follows: 


% = 
°T = 

a = 


meter angle (output) between vane and boom axis 

total angle of attack between resultant air flow and vane 
axis including effects of bending and yaw motion 

angle of attack of rigid missile 

restoring angle of vane 


^ Y f (XqO = elastic axis angular displacement 


cp = attitude error of rigid vehicle. 


The angle-of-attack meter equation of motion is 


cp + 2t a) 9 - 2f cj a + co 2 a = 0, 
b A era b m a v a a * 


(D5) 
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where 


£ = aerodynamic damping ratio 

= mechanic damping ratio of meter 
= frequency of meter. 


Substituting equations . (Dl) through (D4) into (D5) and dividing 
through by gives 


a 


(j 0 


a 


a 

+ 2(t; + t ) ~ + a = a + 

=»m co a v 


X Y 'i (x a ) ’ll 


2? a Y i (x a ) \ . . V Y ’ ( V .. 


ta 


a 


hi + 


‘•fct 


2“ hi 


(D6) 


cp /2 1. X - X 

w a vfc V 


9 


Choosing the sensor frequency » 1, then 


a = a 

v 



Y* (x ) 
1 a 



(x 


Y . (x ) T1. + 

i or ‘i 


£L 


- x ) 

a 


9 


Op. 


(D7) 
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APPENDIX E 

ACTUATOR - ENGINE TRANSFER FUNCTION 

The differential equation relating the output of the actuator to 
the commanded input is given by 


A «o + V 

Kl K„ <f>A + E P A + F »A + G f»A> ' «L A + B L + "l 

where 

P c = actuator command 
P^ = actuator output 
P = control engine gimbal angle 
K q = effective hydraulic spring constant 
K x = open loop gain 
K 2 = valve pressure feedback gain 
= effective load mass 
B = real damping at gimbal 

= effective load spring constant 
A = actuator piston area. 

The transfer function across the actuator is then 

S 3 + ES 2 + FS + G’ 


(El) 


(E2) 


where 


Ki Kq K l 
A(K q + K l ) M l 


(E3) 
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E 


(E4) 


B L , fL K o , K * K o 

-L < K o + K L ) A " A < K o + V 


F 


\ K o 

K + v \ 


+ 


K 3 \ K P \ , ^ K , \ 

«o + V A “ M L A < K o + V «l 


(E5) 


K l K q 1^ 

A < K o +V \ ’ 


(E6) 


The differential equation relating the engine gimbal angle and actuator 
output is 


V + B l P + V =■ K L f3 A 


and has the transfer function 


P 

"3S— 

Pa 


h'*L 


s + 7 s + 

M t 




(E7) 


(E8) 


The detailed block diagram of the actuator is shown in Figure 55. 
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APPENDIX F 


SUMMARY OF VEHICLE DYNAMICS INCLUDING BENDING AND SLOSH 

The equations included in this appendix are written for a vehicle 
in the X, Z coordinate system. The' X axis is chosen tangent to the 
reference flight path moving with the velocity and acceleration of the 
longitudinal vehicle translation. Choosing the coordinate system in this 
manner eliminates the degree of freedom along X whose influence on the 
stability behavior of the vehicle is considered negligible. 

For stability investigations of the vehicle motion, only a short 
interval of flight time is taken. The rotation of the standard flight 
path direction is neglected, and all vehicle parameters (thrust, mass, 
aerodynamic, ect.) are assumed constant. 

The degrees of freedom of the vehicle are 

1. lateral translation of the vehicle rigid body, Z 

2. rotation of the vehicle rigid body relative to an 
inertial system, cp 

3. rotation of the swivel engine relative to an inertial 
system, \|/ 

4. displacement of k propellant slosh masses, Z s ^. 

5. amplitude of i bending modes, r\^ 9 and 

6. the swivel engine compliance, (3 - f3 c . 

The total number of the degrees of freedom for this system is then 

i + k 4* 4 , 

The following equations can be derived by application of Lagrange's 

equation to the kinetic, potential, and dissipative energy of the system 

as was done for the bending and slosh equation in Section V, 
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Angular Acceleration about Center of Gravity . 


cp f = 9 + 9 b + 9 S + 9 e 

9 = " J N a ( x ) a ( x ) dx Cg ' Y 

O 

$3 - T Z{ (X Cg ' X sk> Z S k + (fir) Z 3k} m 3k 

k 


• • 




1 _ 

I 


{ 


[I E + Vl ( 


X 

eg 





<pk 


■T IfE tX cg 


‘ x h + V VV - tI E + Ve <x cg - V 1 Y L (x h ) 



) { Y i <X h> 


(x - X ) Y 1 . (x ) tj . 

v eg h i N h J 'i 


whe r e 


F = F. + F . 
r c 


Acceleration Normal to Vehicle Reference at Center of Gravity 


z' = Z + z\ 4- z 
f b s 


Z 



m 


£ 


Cl ( x ) a( x ) dx 
iN a 
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2 b ■ ■ t I [Y i ( v + "e n<vi\+ 1 1 w\ 


*. - - 10 ?) h, 

k 


Z E V M 


P. 


Acceleration Normal to Vehicle Longitudinal Axis at Sensor Location, 


x 

a 


:<* > = Z - £ J> + 


X W^l 


l*iW 

i 


Bending Equation for Bending Mode 

ri. + 2£. uj. rj. + u? T ). = i- ( F + F + F + F + F ) 
b i x i i M. r\ Z cp a (3 


F 

-a = 

M. 

l 


JT I { m E W 'W - “ E W + *E W Y i (X h> 

1 j 


+ ve n< x h> v j ( y}v ¥ I {w y ;<v + w yv 


j 


•h w y ; 
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(w z sk- s n< x s > z sk 

k 


57 - -^{w- Y i ( v} z 


fk = -L 

M. M. 
1 1 



eg 


+ i E > vv - Weg Y i (x h> - y e w} $ 


+ s 7 { Y i (x h ) - h Y L (x h ) }'i > 


Xi 

et = m" / V x> “ <x) Y i (x) dx 


^ ■ it Ive Y i (x h> - y e Y i (x h ) } e + HT { R,Y i <x h> - Ve * Y I (x h>j- p 


1 1 


where 


g = 


F - D 
M 


M i = / m f (x) [Y i (x )] 2 dx (modal mass), 
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f> f * P+£\i W 
t = 9 + P f 

q: = cp + ~ 
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V 

% = — (yaw plane) 


V cos X 

a w = tan_1 V W - V Sinx (pitch plane) - 

w c 
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APPENDIX G 
BENDING MOMENT 


The. control of large boosters to maintain or reduce the bending 
moment to acceptable levels is frequently an important consideration in 
the design of control systems. The bending moment is produced as a 
result of both external forces and inertial forces acting on the vehicle. 
The structural members of the vehicle must be designed to carry this 
bending moment a-s well as to maintain a structural margin of safety. 

In deriving the expression for the bending moment, all forces are 
resolved parallel and perpendicular to the vehicle longitudinal axis. 
Since the parallel force components act along the longitudinal axis, 
they contribute nothing to the bending moment. The bending moment con- 
tribution from engine dynamics and slosh dynamics is neglected, and a 
rigid vehicle is assumed. 

The free body diagram with force components resolved normal to the 
vehicle axis is shown below. M(x^) represents the bending moment at 
station x^ carried by the structural members. 



FIGURE 59. FREE BODY DIAGRAM OF VEHICLE IN YAW PLANE 
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The forces acting on an element dx are the aerodynamic normal 
force , 


dN = qs [C^ (x)] a dx (Gl) 

and the inertial force 

dP = m f (x) z(x) dx. (G2) 

The prime in the two equations above denotes a derivative with respect 
to x. 


The moment at station x^ due to dN, dP and R 1 is 


x.-x 

l b 


M(x) = R , p(x i - x h ) + 


[dN + dP] (x^ - x) . 


(G3) 


After expressions (Gl) and G(2) are substituted into (G3) , 


M(x) 


R' P(x ± 


x h> + 



[cl! (x)]a + m 1 (x) 
iN a 





x) dx . 


The acceleration z(x) in equation (G4) can be written in terms of OC and 3, 
s ince 


z(x) = z - (x - x) cp: 
eg eg 


(G5) 


and, from equations (39) and (40), 




(G 6) 


(x - x )N' (x - x, )R' 

q> = C S T C P - a - c & y ■ ^ (3. 


(G7) 
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Therefore , 


(x - x) (x - x )N' 1 r n , (x - x) (x - x,)R' 

S(x) - Sl + -£8 i_£S a+ + — £S T C S ^ 

v J m I m I 


Now writing equation (G4) in terms of a and (3 using equation (G8) gives 


M(x) = R' p(x - x) + 


— m' (x) (x - x) dx 

m 1 


(x - x )N ! 


+ — ££ ££ m f ( x ) (x^ - x) dx + qs (x) - x) dxj- a 


X. -X, 
i b 


pi (x - x )R' 1 

— m 1 (x) (x . - x) dx + — m f (x) (x . - x) (x - x) dx f 

m i I i eg J 


Factoring the quantities outside the integral in equation (G9) which 
are independent of x results in a bending moment equation of the form 


M(x ) = M (x )p + M (x )a, 

l pl OC i 


(G10) 


where 


x . -x, 
l b 


M (x . ) = R' (x. - x, ) + — 
6 i i fr m 


m r (x) (x^ - x) dx 


<X M - X h )R ’ 


V*b 


m* (x) (x^ - x) (x - x) dx. 


(Gil) 
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M Q! (x i ) = T J m ' (x) (x i ■ X) dx + qs J 

*b *b 


C* <x) (x - x) dx 
% 1 


+ 



x )N* 
_£E 


x i" x b 

J m' (x) (x t - x) (x cg 
X b 


x) dx* 


(G12) 


Given the aerodynamic characteristics, the mass distribution, the 
vehicle geometry and the control thrust per unit rotation (R ! ), the coef 
ficients Mq, and Mp can be evaluated over the length of the vehicle at 
discrete times of flight* 


Typical values of Mqj(x) and Mp(x) are shown below versus the sta- 
tion location e 


Station In met:«ra , x 



FIGURE 60. Mq,(x) AND Mp(x) FOR A SATURN CLASS VEHICLE AT t = 93 SECONDS 



When a bending moment expression is required in making a control 
study (i.e., shaping the ot and f3 response to reduce bending moment), the 
point of maximum bending is the one of interest. For the vehicle above, 
this would occur in a region around station 22 meters depending on the 
time history response of a and (3. Assuming the response was such that 
the bending moment is maximum at station 22 meters, the expression in 
this case would be 

K - m 

M = - (3. 17(3 + .57a) x 10 7 — . (G13) 
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